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Abstract. Greisen & Calabretta (2002) describe a generalized method for specifying the coordinates of FITS data samples. 
Following that general method, Calabretta & Greisen (2002) describe detailed conventions for defining celestial coordinates 
as they are projected onto a two-dimensional plane. The present paper extends the discussion to the spectral coordinates of 
wavelength, frequency, and velocity. World coordinate functions are defined for spectral axes sampled linearly in wavelength, 
frequency, or velocity, linearly in the logarithm of wavelength or frequency, as projected by ideal dispersing elements, and as 
specified by a lookup table. 
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1. Introduction 

The present paper is the third in a series of papers that estab- 
lishes methods by which information about the physical coor- 
dinates of FITS data may be transferred along with the binary 
image, random groups, and table data. In "Paper I" Greisen 
& Calabretta (2002) describe a revised method for transferring 
coordinate information in the FITS header and outline some 
rules governing the values assigned to the new standard header 
keywords. In "Paper II" Calabretta & Greisen (2002) specify 
the conventions necessary to define celestial coordinates in a 
two-dimensional projection of the sky. This paper defines the 
parameters and conventions needed to specify spectral infor- 
mation including frequency, wavelength, and velocity. In addi- 
tion to these basic conversions, a world coordinate function is 
defined to describe ideal optical dispersers of several common 
types. Several concepts that were suggested by spectroscopic 
issues are generalized to apply to all types of coordinate axes. 
These are a generalized description of non-linear algorithms, 
the -LOG and -TAB non-linear algorithms, and an axis-naming 
keyword. In "Paper IV", Calabretta et al. (2005) specify meth- 
ods to describe the distortions inherent in the image coordinate 
systems of real astronomical data. 

Paper I describes the computation of the world or physical 
coordinates as a multi-step process. The vector of pixel offsets 
from the reference point is multiplied by a linear transforma- 
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tion matrix and then scaled to physical units. Mathematically, 
this is given by 

N 

X j = Si q t = Si m) (Pj ~ , (1) 
j=i 

where pj are pixel coordinates, r, are pixel coordinates of the 
reference point given by CRPIX j, is a linear transformation 
matrix given either by PCLj* or CDLj", N is the dimensionality 
of the WCS representation given by NAXIS or WCSAXES, and s, 
is a scaling given either by CDELT/ or by 1.0. 

The final step in the computation is the conversion of these 
linear relative coordinates into the actual physical coordinates. 
The conventions to be applied to spectral axes, i.e. to frequency, 
wavelength, and velocity axes, are described in this paper. 

2. Coordinate definition 

At this stage in the discussion, we consider the spectral world 
coordinate at the reference point on all other axes. This covers 
the relatively simple case in which the spectral world coordi- 
nate is not dependent on the other world coordinates. Methods 
to describe deviations from this assumption due to the choice 
of spectral reference system are discussed in Sect. 7, while in- 
strumental distortions are discussed briefly in Sect. 9 and, in 
more detail, in Paper IV. 

Spectral coordinates are commonly given in units of fre- 
quency, wavelength, velocity, and other parameters propor- 
tional to these three. The coordinate types discussed here are 
then frequency, wavelength, and "apparent radial velocity" de- 
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True radial velocity / c 

Fig. 1. Conventional velocities as a function of true radial velocity for selected values of the transverse velocity. The apparent radial velocity 
v/c is plotted in the left panel, the radio velocity V/c is plotted in the center panel, and the optical velocity Z/c is plotted in the right panel. 
Note that, at a transverse velocity of 0.9c, the red shift z always exceeds 1. Note that each family of curves intersects the axis of zero velocity 
measure at the same values of v t because the redshift is zero for these combinations of v t and u r . Thus, regardless of the value of v t , for any 
given redshift the velocity measures agree on whether the object appears to be receding or approaching. However, the object may appear to 
be receding at high speed even when it is actually approaching at high speed, though the converse is never true. At transverse velocities of 
vjc > V2/2(« 0.71) all of the velocity measures are positive (receding) regardless of whether the object is actually receding or approaching. 
Furthermore, for vjc below -0.5, the faster the object approaches, the smaller the transverse velocity required to make it appear to be receding! 



noted by the symbols v, A, and v. There are also three con- 
ventional velocities frequently used in astronomy. These are 
the so-called "radio" velocity, "optical" velocity, and redshift, 
denoted here by V, Z, and z and given by V = c(v ( ) - v)/Vo, 
Z = c(A - Ao)/Aq and z = Z/c. The velocities are defined so that 
an object receding from the observer has a positive velocity. 
We assume throughout that each image has at most one spec- 
tral coordinate axis; a similar limit on celestial coordinates was 
assumed in Paper II. 

As discussed by Lindegren and Dravins (2003), the appar- 
ent radial velocity to be described here is itself a conventional 
velocity. The shift of a spectral line from its rest frequency is 
caused by a variety of effects, not just the Doppler shift due 
to a radial motion. The time dilation of an object moving with 
respect to the observer causes a spectral shift, even if that mo- 
tion is entirely transverse. Gravitational fields at the radiating 
object and along the line of sight to the observer will shift the 
observed frequency. Furthermore, the apparent spectroscopic 
velocity may be produced at a peculiar point in the object, 
e.g. upwelling convective cells, rather than at the center of mass 
of the object. It may also be affected by absorption in interven- 
ing material and in particular by the cosmic redshift. We will 
use the term "apparent radial velocity" here simply to refer to 
that pseudo-physical velocity described by the equations pre- 
sented here. The apparent radial velocity then is a sum of all 



spectroscopic effects presented as if they were solely a radial 
velocity. While this may be sufficiently accurate for many uses, 
observers wishing to achieve very high accuracies should con- 
sult Lindegren and Dravins, and references therein, closely. 

The effect of a transverse velocity deserves a little more 
discussion. As given by Lang (1974), 

A = ^o . C + Vr . (2) 

where v x is the true radial velocity, v t is the true transverse ve- 
locity, c is the speed of light, and miscellaneous effects such 
as gravitational redshifts are assumed negligible. Even at or- 
dinary astronomical velocities, this has a measurable effect; a 
transverse velocity of 100 km s _1 produces an apparent velocity 
of 15 ms -1 . The effects of transverse velocities such as might 
be found in astronomical jet sources on the apparent radial, ra- 
dio, and optical conventional velocities are illustrated in Fig. 1 . 
When the transverse velocity is zero, radio velocity ranges from 
-oo to +c, optical velocity from -c to +oo, and apparent radial 
velocity is v t when v t = 0. Because of all of the other processes 
that shift the observed frequency, it is inappropriate to propose 
keywords to describe true velocities at this time. Additional dis- 
cussion of this matter is deferred to Appendix A. 

For the purposes of this section and the next, we con- 
sider those spectral axes at a single celestial coordinate that are 
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Table 1. Spectral coordinate type codes, (characters 1-4 of CTYPE&a). 



Code 


Name 


Symbol 


Associate 
variable 


Default 
units 


FREQ 


Frequency 


V 


V 


Hz 


ENER 


Energy 


E 


V 


J 


WAVN 


Wavenumber 


K 


V 


irr 1 


VRAD 


Radio velocity 


V 


V 


m s' 1 


WAVE 


Vacuum wavelength 


A 


A 


m 


VOPT 


Optical velocity 


Z 


A 


i 

m s 


ZOPT 


Redshift 


z 


A 




AWAV 


Air wavelength 


A* 


A, 


m 


VELO 


Apparent radial velocity 


V 


V 


m s -1 


BETA 


Beta factor (y/c) 


P 


V 





linearly sampled in wavelength, frequency, or apparent radial 
velocity. The radio and optical "velocities" are directly pro- 
portional to frequency and wavelength, respectively, and thus 
do not constitute additional cases for the type of sampling. 
Frequency and wavelength axes may also be regularly sampled 
in their logarithm. Wavelengths are sometimes given "in air" 
rather than in vacuum and denoted here by A- d . This non-linear 
distinction is discussed in Sect. 4. Frequency may also be de- 
scribed in energy (= hv) units and in Kaysers ("wave number" 
= 1/A) units. 

Following the convention of Papers I and II, the first four 
characters of CTYPEfea specify 1 the coordinate type, the fifth 
character is ' - ' and the next three characters specify a prede- 
fined algorithm for computing the world coordinates from in- 
termediate physical coordinates. When k is the spectral axis, 
the first four characters shall be one of the codes shown in 
Table 1 . The mathematical symbols we will use for the codes 
and their default units are also listed. Units are specified with 
the CUNITfea keyword. Standard values for unit strings are de- 
fined in Paper I. The IAU-standard prefixes for scaling the units 
are described in Paper I and should be used with all coordinate 
types, except that the dimensionless ones are not scaled. 

The non-linear algorithm in use is specified by the final 3 
characters of CTYPEfca. In spectral codes, the first of the three 
characters specifies the physical parameter type in which the 
data are regularly sampled, the second character is 2, and the 
third character specifies the physical parameter type in which 
the coordinate is expressed. Non-linear algorithm codes, the 
last four to be introduced later, are summarized in Table 2. 
When the algorithm in use is linear, the last four characters 
of CTYPEfea shall be blank. The original FITS paper by Wells 
et al. (1981) contained a number of suggestions for the values 
of CTYPE ;. Those suggestions are to be superseded by the con- 
ventions of this paper and of Papers I and II. 

The relationships between the basic physical quantities v, A, 
and v are shown in Table 3. The symbols Ao and v are the 

1 While Papers I and II use the generic intermediate axis number 
we use here the axis number k as the spectral intermediate axis num- 
ber. The single-character version code a was introduced in Paper I. It 
has values blank and A through Z. 



Table 2. Non-linear algorithm codes, (characters 6-8 of CTYPEfca). 





Regularly 




Code 


sampled in 


Expressed as 


F2W 


Frequency 


Wavelength 


F2V 


Frequency 


Apparent radial velocity 


F2A 


Frequency 


Air wavelength 


W2F 


Wavelength 


Frequency 


W2V 


Wavelength 


Apparent radial velocity 


W2A 


Wavelength 


Air wavelength 


V2F 


Apparent radial velocity 


Frequency 


V2W 


Apparent radial velocity 


Wavelength 


V2A 


Apparent radial velocity 


Air wavelength 


A2F 


Air wavelength 


Frequency 


A2W 


Air wavelength 


Wavelength 


A2V 


Air wavelength 


Apparent radial velocity 


LOG 


logarithm 


Any 4-letter coordinate type 


GRI 


detector 


Any from Table 1 


GRA 


detector 


Any from Table 1 


TAB 


not regular 


Any 4-letter coordinate type 



rest wavelength and frequency, respectively, of the spectral line 
used to associate velocity with observed wavelength and fre- 
quency. The relationships between the derived quantities and 
the basic quantities with which they are associated are shown 
in Table 4. Derivatives are included in both tables since they 
will be needed in the coordinate computation. 

3. Basic coordinate computation 

This section describes the computation of coordinates that are 
sampled linearly or logarithmically in a coordinate of similar 
type, and of coordinates that are sampled linearly in a coordi- 
nate of different type. The intermediate world coordinate for 
spectral axis k, given by Eq. (1), will be denoted, for conve- 
nience, by 

w = x k . (3) 

The final world coordinate will be denoted by S . At the refer- 
ence point, it has the value 5 r , which is defined by the keyword 
CRVALfea. 

3. 1 . Linear coordinates 

Linear coordinates are represented in CTYPE ka by characters 1- 

4 containing any code in column one of Table 1 and characters 
5-8 blank. The world coordinate value for spectral axis k is 
computed with the above definitions and Eq. (1) as 

5 = S T + w. (4) 

As a general rule, non-linear coordinate systems will be 
constructed so that they satisfy Eq. (4) to first order at the ref- 
erence point. 
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Table 3. Basic spectral transformation equations. 



c 



(6) ^ = _£. 

cU A 2 



V = V ' 



Vc 2 - V 2 



A = A Q - 



Vc 2 ^ 2 



/o-v 2 



(8) 

(10) 
(12) 

(14) 
(16) 



dv 
di; 

dA 

d« 

du 
dv 

do 
dA 



cv 



(c + v) Vc 2 - V 2 

c 
'v 2 



(c - v) Vc 2 - V 2 

Acw\ 
"(v 2 + vg) 2 

4cAA 2 
U 2 + A 2 ) ) 2 



(7) 
(9) 

(11) 
(13) 

(15) 

(17) 



3.2. Logarithmic coordinates 

Data are often sampled in logarithmic increments in one or 
more coordinates. For example, spectra are sometimes sampled 
in logarithmic increments of wavelength or frequency. With 
this type of sampling, the source motion, expressed as V, Z, 
or z, shifts the "pixel" coordinates of spectral features by the 
same amount over the whole image. This allows relative veloc- 
ities between spectra to be determined using cross-correlation 
methods in the pixel arrays. 

For logarithmic axes, the last four characters shall be 
'-LOG'. While there are only three logarithmic coordinates 
commonly used in spectroscopy, namely FREQ- LOG, WAVE - LOG 
and AWAV-LOG, it would be unwise to forbid any coordinate 
type with the -LOG non-linear algorithm. Many such combina- 
tions may have little physical meaning or be intractable math- 
ematically, but these are simply reasons to be cautious when 
using -LOG. This algorithm is evaluated simply 



S v e" 



(5) 



This form of the logarithm satisfies the requirement that 
g^| r = 1 so that Eq. (4) is satisfied to first order at points near 
the reference point. Thus S ~ S r + w. This form of the logarith- 
mic algorithm also has the desirable attribute that the units of 
the coordinate (S), reference coordinate, and scale are the same 
and are of a simple form. The units for CRVALfea, CDELTfea, and 
CDk_ja are specified by the CUNITfca keyword as defined in 
Paper I. These quantities are in normal, non-logarithmic units 
such as 'Hz ' for FREQ-LOG and 'm' for WAVE-LOG. The pre- 
fixes and alternate units described in Paper I may be used, such 
as 'MHz' and 'Angstrom'. 

Logarithmic quantities are frequently expressed as log 
base 10 rather than as natural logarithms. For such data, the 
CDELTfea and CDk_ja will need to compensate by including a 
factor ofln(10). 

3.3. Coordinate axis names 

The generality of this algorithm, and the -TAB algorithm to be 
introduced below, suggest the need for a more general descrip- 



Table 4. Extended spectral transformation equations. 



v = v (l ) 

c 



dv vn 

(18) AW=-- 

dV c 



V = c 



vq - v 
v 



v = E/h 
E = hv 

V = CK 
K = v/f 



A = A (l + -) 

c 



(20) 
(22) 
(24) 
(26) 
(28) 



dV 
"dV 

dv 
dE 
dE 
dv 
dv 

dK 
dK 

dv 



c 



— = l/h 



— = l/c 



(30) — = - 



dA _ Ao 
dZ ~ c 



Z = c 



A-A Q 



A 

A = A (l+z} 
A-A 



(32) 



dZ _ c 
dA~ 1^ 

(34) ^=A 



(36) 



dz 

dz 
dA 



1 



(38) — = 



dv 
dp 



P = vie 



(40) 



dp 
dv 



l/c 



(19) 
(21) 
(23) 
(25) 
(27) 
(29) 

(31) 
(33) 
(35) 
(37) 

(39) 
(41) 



tion of the coordinate than may be indicated in the first four 
letters in the value of CTYPEfea. We hereby reserve the keyword 



CNAME ia 



(character-valued) 



to provide a description of a particular coordinate in up to 68 
characters. Its default value will be all blank. For binary table 
vectors, the keyword will be iCNhna, while for pixel lists it will 
be TCNAna, where ;' is the intermediate world coordinate axis 
number and n is the table column number to which the key- 
word applies. This keyword may be used with standard axis 
types such as GL0N/GLAT or FREQ, but will be of greatest value 
with non-standard axis types such as, hypothetically, CTYPE i = 
'TFPY-LOG' meant to indicate "log of Y position in telescope 
focal plane," which would be recorded in the CNAME ia card for 
that axis. This keyword provides a name for an axis in a partic- 
ular WCS, while the WCSNAMEa keyword names the particular 
WCS as a whole. 



3.4. Non-linear spectral coordinates 

We now consider the case where an axis is linearly sampled in 
spectral variable X, but is to be expressed in terms of variable 
S. 
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3.4.1 . Spectral algorithm codes 



3.4.2. Spectral algorithm chain 



Given the large number of spectral variables in Table 1 it 
is clear that there are very many pairwise combinations; in 
general, the relationship between any pair may be non-linear. 
However, each of the spectral variables in Table 1 is linearly 
related to one or other of v, A, A a , or v. Thus we may restrict X 
to one of these four basic variables. 

Even with this restriction on X there are still many possi- 
ble combinations of X and S . In order to reduce the number 
still further we introduce an intermediate variable, denoted by 
P, that is the basic variable, v, A, A. d , or v, with which S is as- 
sociated via a linear relation. This associate variable is listed 
for each spectral variable in column 4 of Table 1 . Thus the se- 
quence of transformations may be summarized as 



Pj — * Xk (= w) — > X ~» P 



where the non-linear transformation, indicated by the wiggly 
arrow, is between X and P. 

Table 3 lists the equations, X = X(P), and their inverses, 
P = P(X), linking the basic spectral types v,A, and v (dis- 
cussion of A. d is deferred to Sect. 4). These equations are gen- 
erally non-linear. Likewise, Table 4 lists the linear relations, 
S = S(P), and their inverses, P = P(S), linking each spectral 
variable with its associate. When S is one of the basic types, 
v, A, /l a , or v, then P = S and S (P) is identity. 

Thus the functional relationship between S and X is speci- 
fied via intermediate variable P as S (X) = S (P(X)) with inverse 
X(S ) = X(P(S )). Since S (P) is linear, P and X must always dif- 
fer, otherwise P(X) would be identity with the result that S (X) 
would be linear, contrary to the assumption of a non-linear axis. 

Non-linear spectral coordinate codes are constructed on 
this basis by combining the spectral coordinate type code for 
S from Table 1 with the non-linear spectral algorithm code 
in Table 2. The first letter of the algorithm code defines X as 
frequency (F), wavelength (W), air wavelength (A), or appar- 
ent radial velocity (V), and the third letter likewise defines P. 
For example, in Z0PT-F2W, X is Frequency, P is Wavelength, 
and the non-linear conversion between the two (F2W) is de- 
fined in Table 3, Eq. (10). The desired spectral coordinate S is 
redshift (ZOPT), and this is related to associate variable P, the 
Wavelength, by the linear equation given in Table 4, Eq. (36). 

It is apparent from the foregoing that it is possible to con- 
struct invalid spectral CTYPE coordinate codes. For example, 
Z0PT-F2V is unrecognized since z is associated with A, not v. 
That is not to say that this code could not be interpreted in 
principle — after all, equations linking z and v could have been 
included in Table 4 — simply that it is not recognized in prac- 
tice. The associate variable, P, was introduced in the first place 
to reduce the possible number of combinations and adding new 
ones like this would defeat that purpose. This is particularly 
relevant to the software implementation; Z0PT-F2V would tell 
software to chain its v - v function with av-z function but, in 
general, the latter function will not have been defined. 



Consider now the computation of spectral coordinate S for an 
axis that is linearly sampled in X. The statement that an axis is 
linearly sampled in X simply means that 



X = X x + w 



dX 

dw' 



(42) 



where dX/dw is constant. This constant is determined by im- 
posing the requirement that 



dS 
dw 



1, 



(43) 



so that Eq. (4) is satisfied to first order at points near the refer- 
ence point: 



S ~V 

Thus 
dX 
dw 



w. 



dP 
dS 



dP 

1 dX 



(44) 



(45) 



which gives dX/dw as a function of X r , dP/dS being constant 
and dP/ dX a function of X. Given that the functions S — S (P) 
and P = P(X) are known, as are their inverses X = X(P) and 
P = P(S), then the equation for S as a function of w may be 
obtained from Eqs. (42) and (45) 



S(w) =S\P\X(P(S T ))+w 



dP 
dS 



I 



dP 
dX 



(46) 



where S T is given by CRVALfca. 

Equation (46) suggests a three-step algorithm chain: 

1. Compute X at w using Eq. (42); X r = X(P(5 r )) and dX/dw 
are constants that need be computed once only, the latter 
obtained from Eq. (45) as a function of X r . 

2. Compute P from X using the appropriate equation from 
Table 3. 

3. Compute S from P using the appropriate equation from 
Table 4. 

The inverse computation by which the intermediate coordinate 
w is computed for a given value of S is effected by traversing 
the algorithm chain in the reverse direction. The inverse equa- 
tions required are all listed in Tables 3 and 4. The inverse of 
Eq. (42) is, trivially, 

dX 

w = (X-X t )/~ (47) 
dw 

3.4.3. Example non-linear calculation 

As an example of a non-linear coordinate computation, con- 
sider a Z0PT-F2W axis. F in the F2W code indicates that the axis 
is linearly sampled in frequency, but it is to be expressed in 
terms of redshift as indicated by the spectral coordinate type of 
ZOPT. Table 1 indicates that redshift is associated with wave- 
length, hence the W in F2W is correct. Of the X, P, and S vari- 
ables in the preceding section, the axis is linear in frequency so 
X is v, the associated variable is wavelength, so P is A, and the 
required variable is redshift, so S is z. 
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Table 5. Sample non-linear coordinate combinations, where w, the intermediate coordinate of Eq. (3), has a different meaning for each of these 
equations It has units given by the spectral coordinate type code as listed in Table 1. 



FREQ-W2F v = 



FREQ-V2F y : 



v r - w 



Vr 



VEL0-W2V v = c 



VRAD-W2F V = 



WAVE-F2W /I : 



V r (c ~ Vr) + CW 

c-V r + w 

A, 



VEL0-F2V v = c 



A r — W 

Al-l 



V0PT-F2W Z = 



Z0PT-F2W z ■- 



A 2 + l 

XAc + Z r ) + cw 
c + Z r - w 

Zr(l + Z r ) + W 
1 + Zr - W 



WAVE-V2W A = An 



V0PT-V2W Z = , 



1-S, 



1 + S V 



c 2 - v 2 + CW 

(C - V r ) -y/c 2 " v f 



(48) 

(49) where A,, 

(51) 

(52) 

(53) where A v 

(55) 

(56) 



(57) where B A = 2 ° ° 



(c + v r ) -^c 2 - V 2 



(59) where fi„ 



1 (61) where £» z = 



(Ai + A?y 

(vg + v 2 ) 2 

y| - 1 + 4y z w/c 

(1 + ¥ 2 Z ) 2 



(50) 



(54) 



(58) 
(60) 

(62) with y z =i+z r /c 



(63) 



Since CRVALfea for the Z0PT-F2W axis would be recorded 
as a redshift, z r , this must first be converted to frequency by 
applying Eqs. (34), 

A r = /to(l+Zr) 

and (6), 

c c 

A x Ao(l+Zr) 

This provides X r for Eq. (42). Then dv/dw (i.e. dX/dw) is ob- 
tained from Eqs. (35) & (1 1) evaluated at the reference point: 



dv 
dw 



dA 
di" 



-c v 2 

Vf VQ 



Redshift may now be computed for any given value of w. The 
first step is to compute v (i.e. X) at w using Eq. (42): 



V r + W ■ 



dv 
dw 



In this instance w will be a redshift. Then A (i.e. P) may be 
computed from v using Eq. (10) from Table 3: 

A = c/v. 

The third and final step is to compute z (i.e. S ) from A using 
Eq. (36) from Table 4: 



z = 



^0 



It is also possible to combine the three steps into a single 
equation, although there are many pair-wise combinations of 
spectral variables and for some the combined equations may 
be quite cumbersome. For the sake of illustration, consider the 
above equations for the Z0PT-F2W axis. Eq. (42) becomes 

c cw 

V ~ A (l+Zr) ~ ioU+Zr) 2 ' 

Substituting this in the equations for A and z and simplifying 
we obtain 



z = 



Z r (l +Z r ) + W 
1 + Zr + W 



A representative sample of these direct translation equations is 
shown in Table 5. They are useful in showing the nature of the 
coordinate non-linearity. 

One thing to notice about the equations of Table 5 is that 
the meaning of w differs for each. For example, Eq. (55) may 
not be obtained from Eq. (56) simply by multiplying both sides 
by c: in Eq. (56) w is a redshift (ZOPT), whereas in Eq. (55) it 
is an optical velocity (VOPT). 

3.4.4. Coordinate parameters 

Aside from CRVALfca, the coordinate computations of Sect. 3.4 
require one extra parameter when evaluating the F2V, V2F, W2V, 
V2W, A2V, and V2A non-linear algorithms, namely the rest fre- 
quency or wavelength of the spectral-feature of interest. These 



are fundamental physical parameters so, rather than use the 
PVijna parameters defined in Paper I, which would tend to dis- 
guise them, the special floating-valued keywords 



and 



RESTFRQa (floating-valued), 



RESTWAVa (floating-valued), 



are hereby reserved for the purpose. They are represented by 
symbols vo and Aq, respectively. Their units are ' Hz ' and ' m ' 
respectively, fixed to save having additional keywords to define 
them. RESTWAVa is for wavelengths in vacuum only. One or 
the other of these keywords must be included for the above- 
mentioned algorithm codes; usually RESTFRQa would be given 
for F2V and V2F, and RESTWAVa for the others. FITS writers 
should always write one or other of these when it is meaningful 
to do so, even for algorithm codes such as F2W or W2A that do 
not require them. 

Keyword RESTFREQ has been used in previous FITS files 
and should be recognized as equivalent to RESTFRQ. 

4. Air wavelengths 

The wavelengths so far discussed are measured in vacuum. 
However, dispersion coordinates for UV, optical, and IR spectra 
at A > 200 nm are commonly given as wavelengths in air. The 
relative difference between the two varies between 0.028% and 
0.032% across this range. To identify wavelengths measured 
in dry air at standard temperature and pressure rather than in 
vacuo, we introduce the coordinate type AWAV for which the 
reference value and increment must be expressed accordingly. 
The two measures of wavelength are related by 



A = n(A & )A a , 



(64) 



where A is the wavelength in vacuum, A a is the wavelength in 
air, and n(A a ) is the index of refraction of dry air at standard 
temperature and pressure. This varies non-linearly with wave- 
length. The standard relation given by Cox (2000) is mathemat- 
ically intractable and somewhat dated. The International Union 
of Geodesy and Geophysics (1999) adopted the rather simpler 
formula 2 

„(^= 1 + 10- ( 287.6155+ ™ + ^), (65) 

where A d is the wavelength in micrometers. This formula suf- 
fices for conversions from air to vacuum when no more than 
0.25 parts per million accuracy is required. The derivative, 
which may be required in Eq. (45), is 



dA_ 

di a 



= 1 + lO" 6 287.6155 - 



1.62887 0.04080 \ 



(66) 



While the inversion of Eq. (65) is algebraically intractable, 
the vacuum wavelength may be used to evaluate Eq. (65) since 

2 The quoted formulas apply only at normal optical wavelengths. In 
the UV and IR spectral domains, atmospheric absorption lines cause 
refractivity to be a strong, weather dependent function of frequency. 
See, for example, Mathar (2004). 



it differs so little from the air wavelength. The resulting error in 
the index of refraction amounts to 1 : 10 9 , and this is less than 
than the accuracy of the empirical formula. Thus, 



A/n(A). 



(67) 



As usual, an axis that is linearly sampled and expressed in 
air wavelengths is described with a CTYPEfca of 'AWAV and 
evaluated with Eq. (4). 

Algorithm codes for the non-linear conversions are A2F, 
A2W, A2V, and their complements, F2A, W2A, and V2A as 
listed in Table 2. Use of the three-step procedure described in 
Sect. 3.4 would require A. d as a function of each of the other 
spectral variables, together with their inverses and derivatives. 
Thus it is much simpler to handle air wavelengths as a sepa- 
rate, extra step in the algorithm chain. For example, to compute 
world coordinates for VRAD-A2V, the value of CRVALfea would 
first have to be converted from radio velocity to air wavelengths 
via Eqs. (18), (10), and (67), and dAJdw forEq. (45) would be 
obtained from Eqs. (19), (7), and (66). Then, the value of A- d 
computed for w via Eq. (42) would be transformed to vacuum 
wavelengths via Eq. (64), and then to radio velocity via Eqs. (6) 
and (20). 

5. Dispersed spectra 

One common form of spectral data is produced by imaging the 
light from a disperser. The wavelengths of the light at some 
position in the image are related to the position and wavelength 
of the light in the field of view illuminating the disperser; the 
relation is defined by the physics of the disperser. In general 
the light received at a pixel in the detector is a superposition 
of different wavelengths from different points in the field of 
view. However, if the field of view is limited spatially, usually 
by aperture masks and fiber optics, each pixel receives light 
from only a small range of wavelengths. It is then possible to 
define a world coordinate function relating pixel position and 
effective wavelength. This is the basis of many astronomical 
spectrographs. 

In the following section we use the physical relation appli- 
cable to the dispersers commonly used in astronomical spec- 
trographs to define a world coordinate function for comput- 
ing spectral coordinates. The relation applies to the simple, 
though common, case of single dispersers. More complex spec- 
trographs with multiple dispersers, such as those using multiple 
passes through prisms, are not described by the methods of this 
section. The equations developed below also assume that the 
radiation enters perpendicular to the face of the prism, a condi- 
tion not met by some widely used spectrometers. Alternatives 
to using this ideal world coordinate function, based on the 
physics of simple dispersers, are the table lookup described in 
Sect. 6 and empirical function fits provided by the the distortion 
function mechanism described in Paper IV. 

We require that the dispersion occurs along just one direc- 
tion on the detector and that the intermediate coordinates are 
computed so that only one world coordinate axis corresponds 
to wavelength. Paper IV describes how distortions and effects 
of the aperture shape can be removed to satisfy this require- 
ment. The distortion correction is also used to remove aberra- 
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Fig. 2. Geometry of gratings, prisms, and grisms. This simplified representation omits the collimation and focusing optics. Dashed lines mark 
ray paths in the plane of the figure - the "dispersion plane". The normal to the grating/exit prism face and the normal to the detector plane are 
each projected onto the dispersion plane, and angles a, y, and 9 are measured with respect to these projected normals. Usually the incident 
ray for a prism or grism is perpendicular to the entry face so that a is equal to the prism angle, p. Angle y is wavelength-dependent, and 
consequently so is the offset f in the dispersion direction on the detector. The intermediate spectral world coordinate, w, is proportional to 
Reference wavelength A r follows the reference ray defined by y r and illuminates the reference point at w = tj = 0. The normal to the detector 
plane is shown tilted by angle 6 from the reference ray though typically this angle is zero. The grating spacing G { is indicated. 



tions causing departures from the ideal physical behavior of the 
dispersers assumed here. 

The dispersers we consider are gratings, prisms, and 
grisms, which are a combination of a prism with a grating 
within or on the surface of the prism. Gratings may be re- 
flecting or transmissive and may be fabricated with surface 
relief rulings, holographic surface relief patterns, and holo- 
graphic volume phase patterns. In the discussion we use the 
terms grooves, lines, and ruling to refer to the periodic diffrac- 
tion structures that produce the interference. By combining the 
two physical equations for interference and refraction a single 
world coordinate function covers all these cases with appropri- 
ate choice of parameters. 



5.1. The grism coordinate function 

This section defines a world coordinate function using the basic 
laws of refraction and interference. It is beyond the scope of 
this discussion to give the full background for the laws and 
concepts underlying dispersive spectroscopy. Of the many texts 
on the subject, a standard reference is Astronomical Optics by 
Schroeder(1999). 

Spectroscopic dispersers are based on interference and re- 
fraction effects which are naturally described in terms of wave- 
length. Moreover, they are most often used in the regime where 
wavelength is commonly the spectral coordinate. For these rea- 
sons the grism function is defined in terms of wavelength. 
The function assumes that all of the dispersion occurs at one 
place. This is why combinations of dispersers, other than a sin- 
gle grism, are not described by this function. This assumption 
means that for prisms and grisms the function is rigorously cor- 
rect only when the light enters perpendicular to the face of the 
prism and the grating is at the exit face. However, even when 
these conditions are not exactly met the function may still be a 
good approximation with slight adjustments of the parameters. 



5.1 .1 . The grism equation 

The basic physical relationship between the wavelength, A, the 
angle of incidence of the light, a, and the diffracted/refracted 
angle, y, is given by a combination of the grating interference 
equation and Snell's law of refraction: 3 

GmA 

= n(A) sin a + sin y, (68) 

cose 

where G is the grating ruling density, m is the interference or- 
der, and n(A) is the wavelength-dependent index of refraction 
of the prism material. For a pure prism the order m is zero and 
for a reflection or transmission grating n(A) is the unit func- 
tion. The plane containing a and y is the dispersion plane, il- 
lustrated in Fig. 2, and the angles are measured relative to the 
projection in the dispersion plane of the normal to the grating 
or exit prism face. Usually the normal does lie in the disper- 
sion plane, but small values of e, the angle between the normal 
and the dispersion plane, are sometimes used for instrumental 
design reasons. 

By convention, the angles of incidence and diffrac- 
tion/refraction are measured from the normal to the beam and a 
is always measured positive in the anticlockwise direction. The 
sign changes on either side of the normal and this determines 
the sign of y. Thus, for the reflection grating in Fig. 2, a > 
and y < 0, and vice versa for a transmission grating, prism, or 
grism. 

Reflection grating geometry is sometimes defined by the 
angle <f> measured from the camera axis to the collimator axis 
and the tilt t of the bisector relative to the grating normal. If <p 
and t obey the same sign convention as or, then a = (p/2 + t. 

In spectrographs with prisms or grisms the requirement that 
the incident light be normal to the prism entry face means that 
a is equal to the prism angle (p in Fig. 2). Even for oblique 
entry, identifying a with the prism angle is often a good initial 
approximation. 

3 Also known as Descartes' Law of Sines; see for example, 
http : //wikipedia . org/wiki/Snells_law 
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The requirement that the difFraction and refraction occur at 
one point means that a is fixed and independent of wavelength. 
The variation of y with wavelength then defines the relation be- 
tween wavelength and position £ on the detector. The reference 
wavelength A r is the wavelength at the reference point corre- 
sponding to the zero point of £ and the intermediate spectral 
world coordinate w. 

The prism's dispersive power derives from the variation of 
its index of refraction with wavelength. While this variation 
depends on the material and is generally non-linear, it is suffi- 
cient to approximate it by a first-order Taylor expansion about 
the reference wavelength, A r : 



n(A) « n T + n' T (A - A r ), 



(69) 



where we have written n r as a shorthand for n(A t ), and likewise 
n[ for dn/d/i| r . Combining Eqs. (68) and (69) yields the grism 
equation 

(n r - n' r A r ) sin a + sin y 



A = 



Gm/ cos e-n' T sin a 



(70) 



where the denominator must not be zero, though Gm, n[, and a 
may be zero in different types of spectrographs. 

In order to define a world coordinate function we need A 
as a function of the intermediate world coordinate, w, which 
is proportional to Since Eq. (70) gives A as a function of y 
it remains to determine the relationship between y and First 
note that Eq. (70) evaluated at A = A r provides 



y x = sin 1 (GmA r / cos e - n r sin a) . 



the exit angle of the reference ray for which w = £, — 0, as in 
Fig. 2. 

Fig. 2 shows angle 9, which is measured from the refer- 
ence ray to the camera axis using the same sign convention as 
y, i.e. if y is clockwise-positive as for a grism then so is 8. 
Normally 8 would be zero, but it is included here to provide a 
small correction. Depending on the sign convention for y, sim- 
ple geometry for a flat detector gives 



y = Tr + + tan _1 (±f// - tan 9), 



(72) 



where / is the effective focal length of the camera. The plus 
sign is taken when the sign convention for y is such that 
£ increases as y increases, and the minus sign otherwise. 
Section 5.1.3 discusses the resolution of this potential sign am- 
biguity. 

5.1 .2. GRI coordinate axes 

In keeping with the preceding sections we wish to define gen- 
eral grating, prism, and grism world coordinate representations 
such as WAVE-GRI, FREQ-GRI, etc. 

Bearing in mind that the grism equation, Eq. (70), is ex- 
pressed in terms of wavelength, then given a FREQ-GRI axis, 
for example, it would be straightforward to convert the refer- 
ence frequency, v r , given by CRVALfe, from frequency to wave- 
length via A = c/v. However, it would not be valid to convert 
the intermediate spectral world coordinate, w, from frequency 
to wavelength like this because it is not a true frequency. While 



it may be a close approximation near the reference point it de- 
viates at points away from it. 

Thus, interpretation of an axis such as FREQ-GRI neces- 
sarily involves a procedure similar to that of Sect. 3.4, and to 
apply that methodology we need a parameter, the grism param- 
eter, T, that is a known function of the spectral variables and for 
which the axis is linearly sampled; this will substitute for X in 
Eq. (42). Since £ is proportional to w 

€ = crw, (73) 
where cr is a constant. Combining this with Eq. (72) we have 



tan(y - y r - 9) — - tan 9 ± wcr/ f. 

Thus the grism parameter may be identified as 

r = tan(y -y t — 9), 

and this satisfies Eq. (42) at the reference point: 
T r = - tan 6. 



(74) 



(75) 



(76) 



The grism parameter has a simple geometrical interpretation in 
Fig. 2; it is the offset on the detector from the point of intersec- 
tion of the camera axis measured in units of the effective focal 
length, /. 

Following Sect. 3.4 (with X replaced by F), we have 



dr 

r r + w— 

aw 



py. where dT/ dw is constant: 



dr 
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dP 


dS 


dw 


" dP 


r dS 


i dw 



(77) 



(78) 



As before, S is the spectral variable in which the axis, and 
hence w, is expressed, and P is the basic variable, v, A, A a , or v, 
with which S is most closely associated. Recognized spectral 
variables and their associates are listed in Table 1 . When S is 
one of the basic types, v, A, A a , or v, as is often the case, then 
P = S whence dP/dS| r = 1. 

As previously, we require that CDELT; or CDLj be set so 

that 



dS_ 

dw 



1 



(79) 



(exactly how this is done is addressed in Sect. 5.1.3). Since we 
only have T directly in terms of y we may use 



(80) 



where dF/d/l| r itself may be deduced from the derivatives of 
Eqs. (75) & (70): 



dr 


dr 


dA 


dP 


r ~ dA 


r dPr 



dr 

dA 



dr 

dy 



I 



dA 
dy 



Gm/ cos e-n[ sin a 



cosy r cos 2 9 
Combining Eqs. (78) to (81) we have 



dr Gm/ cos e-n' x sin a dA 
dw cos y r cos 2 8 dP 



dP 
dS 



(81) 



(82) 



1 \J 

Table 6. Grism parameters, their default values, and required units. 



Keyword Default Units 



PVk.Oa = G 





nr 1 


PVkAa = m 







PVL2a = a 





deg 


PVk3a = n r 


1 




PVkAa = n' r 





m- 1 


PVk.Sa = e 





deg 


pv/t_6a = e 





deg 



where dA/dP\ r and dP/dS | r come from Tables 3 & 4 respec- 
tively. In the common case where S , and hence P, is wavelength 
then dA/dP\ T and dP/dS | r are both unity. 

The foregoing provides everything needed to compute F 
from w, and this forms the first step of a five-step algorithm 
chain for computing S from w. 

1. Compute r at w using Eq. (77). T r and dr/dw are constants 
that need be computed once only; T r from Eq. (76) and 
dT/dw from Eq. (82) using the appropriate equations from 
Tables 3 & 4. 

2. Compute y from T using the inverse of Eq. (75): 

y = tan-\r)+y T + 6. (83) 

3. Compute the wavelength A from y via Eq. (70). 

4. Compute P from A using the appropriate equation from 
Table 3. 

5. Compute S from P using the appropriate equation from 
Table 4. 

If S , and hence P, is wavelength then the last two steps are not 
performed. 

As usual, the inverse computation is effected by traversing 
the algorithm chain in the reverse direction. The inverse equa- 
tions required for the backward steps have already been given 
except that for Eq. (70): 

y = skr'(/!(Gm/ cos e - n' r sin a) - 

(n r -«^ r )sina), (84) 

and also the inverse of Eq. (77) which is trivial. 

Grism parameters required for these calculations are pro- 
vided by the PVlcjna keywords defined in Table 6, with y r given 
by Eq. (71), and A r computed from the coordinate reference 
value, CRVALfc, by the appropriate equations from Tables 3 & 
4. Default values for missing parameters are also defined in the 
table. Generally only the first three parameters will appear for 
gratings and the first five for grisms. 

5.1 .3. Determination of the grism scale 

A question raised above was how CDELT; or CDLj may be set 
by a WCS composer so that dS /dw\ T = 1. 
From Eqs. (1) & (3) 

w = s k q k (85) 



so 



dw 


dw 


dS 


dP 


dA 




dS 


r dP 


r dA 


r dq k 



(86) 



Now, dw/dS | r = 1 by design, dS /dP| r and dP/dA\ T come from 
Tables 4 & 3 respectively, and dA/dq k \ r is the wavelength dis- 
persion (e.g. in nm/pixel) measured at the reference point and 
this is an instrumental parameter. Thus we have everything 
needed to compute s k . When S is wavelength, then all but the 
last factor in Eq. (86) is unity, and s k = dA/dq k \ T . 

It is the correct choice of sign for dA/dq k \j. that resolves the 
sign ambiguity in Eq. (74). 

5.1 .4. GRA: grisms in air 

Thus far we have ignored the distinction between vacuum and 
air wavelengths in the discussion of grism world coordinates. 
In fact, the A variable that appears in the grism equation may be 
either, and in general n(A) in Eq. (68) is the index of refraction 
of the prism material with respect to the surrounding medium, 
either air or vacuum. 

If the dispersion takes place in vacuum, as it does, for ex- 
ample, in a spectrograph on a spacecraft, then the grism equa- 
tion is correct for vacuum wavelengths; if in air, it is correct for 
air wavelengths and A in Sects. 5.1.1, 5.1.2, and 5.1.3 should 
be replaced everywhere by i a . Then as discussed in Sect. 4, 
the computation of quantities such as dA d /dP\ T in Eq. (82) is 
best handled by using the vacuum wavelength, A, as an inter- 
mediary, effectively introducing an extra step into the algorithm 
chain. 

In order to distinguish between grisms operated in air and 
in vacuum we hereby reserve the use of GRI exclusively for 
vacuum operation and introduce GRA for operation in dry air at 
standard temperature and pressure. 

Note that for real spectrographs operated in air at an ob- 
servatory, the actual wavelength system that GRA describes is 
for air at the local conditions. The coordinate value and incre- 
ment at the reference point are normally adjusted to the val- 
ues at standard conditions during calibration. If the accuracy 
of the wavelength measurement requires it, any further correc- 
tion between wavelength at local conditions and wavelength at 
standard conditions may be accomplished via the methods of 
Paper IV. 

5.2. AWAV-GRA examples 

This section illustrates application of the GRA world coordi- 
nate function with three real- world examples of spectral images 
from Kitt Peak National Observatory (KPNO) spectrographs 
(Figs. 3-5). Each spectrum is of an arc calibration lamp that 
produces emission lines of known wavelength. The position of 
each line in the image along the spectral world coordinate axis 
is measured by centroiding on the spectral line profile and iden- 
tified with the known rest wavelength to create a list of pixel 
positions and wavelengths. 

Each example figure shows plots of the positions and wave- 
lengths for a particular spectrograph. The pixel positions are 
plotted along the bottom axis. Rather than plot the wavelengths 
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KPNO COUDE: LONG FOCAL LENGTH 
AWAV (Angstrom) 

6000 5800 5600 5400 5200 5000 4800 
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Pixel 

CTYPE1 = ' AWAV-GRA' / Grating dispersion function 

CUNIT1 = 'Angstrom' / Dispersion units 

CRPIX1 = 1801.7 / [pixel] Reference pixel 

CRVAL1 = 5225.2 / [Angstrom] Reference value 

CDELT1 = -0.4334 / [Angstrom/pixel] Dispersion 

PV1_0 = 3.16E5 / [irf(-l)] Grating density 

PV1_1 = 1 / Diffraction order 

PV1_2 = 13.9 / [deg] Incident angle 



KPNO HYDRA: ECHELLE ORDER 11 

AWAV (Angstrom) 

5250 5200 5150 5100 5050 5000 
H -20 




500 1000 1500 2000 

Pixel 

CTYPE1 = 'AWAV-GRA' / Grating dispersion function 
CUNIT1 = 'Angstrom' / Dispersion units 
CRPIX1 = 944.8 / [pixel] Reference pixel 

CRVAL1 = 5136.8 / [Angstrom] Reference value 
CDELT1 = -0.1287 / [Angstrom/pixel] Dispersion 
PV1_0 = 3.16E5 / [m"(-l)] Grating density 
PV1_1 = 11 / Diffraction order 

PV1_2 = 64.8 / [deg] Incident angle 



Fig. 3. KPNO Coude Feed spectrograph with a long focal length and 
a3KCCD. 

directly, where it would be difficult to see the shape of the 
curve, the difference or correction between the known (air) 
wavelengths and the simple linear world coordinate function 
AWAV (A. d = S r + w) are shown. The corrections are plotted 
along the left axes. In these examples the wavelength units are 
consistently Angstroms. 

The top axis is labeled with simple AWAV linear coordinates. 
The right axis divides the wavelength correction by CDELT1, 
the linear dispersion at the reference pixel. This represents the 
shift, in pixels, of the wavelengths on the detector relative to 
where they would occur in a spectrograph with a linear disper- 
sion relation. 

The departure of the data points from zero indicate the mag- 
nitude of the world coordinate errors that would occur by using 
the simple linear AWAV world coordinate function. The solid 
lines in the figures are the difference between the wavelengths 
produced by the AWAV-GRA grating coordinate function and the 
linear AWAV coordinate function evaluated at the pixel positions 
in the image. 

The usefulness of the ideal grating coordinate function is 
that the curves go through the measured data points with an ap- 
propriate choice of parameters. The parameters are essentially 
those known for the spectrograph and disperser with small ad- 
justments to some of the parameters to produce the best fit to 
the calibration data. (These small adjustments may be viewed 
as corrections for the simplifications in the optical model.) In 
these examples the dispersion function represented by the grat- 
ing world coordinate relation is as good as typically provided 
by empirical polynomial functions. Higher order effects due to 
aberrations are corrected by the distortion correction methods 
defined in Paper IV. 

A close reading of the equations above will reveal that 
the seven grism parameters listed in Table 6 are not inde- 



Fig. 4. KPNO Hydra Fiber Bench Spectrograph using an echelle grat- 
ing in 1 1 th order with a 2K CCD. 

KPNO MARS: VPH GRISM 
AWAV (Angstrom) 

6000 7000 8000 9000 10000 




200 400 600 800 1000 1200 1400 1600 

Pixel 

CTYPE1 = 'AWAV-GRA' / Grating dispersion function 
CUNIT1 = 'Angstrom' / Dispersion units 
CRPIX1 = 719.8 / [pixel] Reference pixel 

CRVAL1 = 7245.2 / [Angstrom] Reference value 
CDELT1 = 2.956 / [Angstrom/pixel] Dispersion 

PV1J8 = 4.50E5 / [m"(-D] Grating density 
PV1_1 = 1 / Diffraction order 

PV1_2 = 27.0 / [deg] Incident angle 

PV1_3 = 1.765 / Reference refraction 

PV1_4 = -1.077E6 / [itTC-l)] Refraction deriv 

Fig. 5. KPNO MARS spectrograph with a 450 lines/mm volume phase 
holographic grism and a 2K CCD. 



pendent. We have chosen these parameters because of their 
physical meaning. However, the independent parameters are 
Gm/ cos(e), n r sin(a-), n[ sin(a-), and 6. It is these combinations 
of parameters which must be used in fitting data. 

Each figure shows the arc line measurements, the coordi- 
nate function curve, and the relevant world coordinate key- 
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words used to produce the curve. Not all of the WCS keywords 
are shown. 

The first example shows the behavior of a long focal length 
spectrograph with a reflection grating used at a low angle of 
incidence. Because of the long focal length the deviation from 
linearity is relatively small though still clearly significant. The 
grating has a density of 316 lines/mm and is used in first order 
to produce a spectrum centered near 5225.2 A with a dispersion 
of 0.43 A/pixel. 

The next example uses a 316 lines/mm echelle reflection 
grating, a grating designed for use at large angles from the 
grating normal, operated in a higher order. It is used to pro- 
duce a spectrum centered near 5136.8 A with a dispersion of 
0.13 A/pixel. The departure from linearity is in the opposite 
sense from the other examples because the echelle grating is 
used with the diffracted angle at a greater angle than the inci- 
dence angle. The higher order results in a fairly large departure 
from a linear WCS. 

The last example illustrates use of a grism. However, this is 
not the more common grism with a ruled grating on the output 
face and the input face oriented normal to the incident beam. 
Instead, the grating is sandwiched in the middle of the prism. 
The prism is oriented with the beam entering and leaving the 
prism at equal angles to the faces resulting in a straight through 
configuration as with a standard grism. 

Another unusual feature of this grism is that it uses a vol- 
ume phase holographic (VPH) grating. While the intensity re- 
sponse is different from a surface relief grating (ruled or holo- 
graphic) the dispersion behavior is the same. 

The full optical equation is complex even in this symmetric 
configuration, but as shown in the figure using the GRA func- 
tion, a very good description of the coordinates can be ob- 
tained. 

The prism has a 27° angle with an index of refraction of 
1 .764 near the reference wavelength. The grating has an equiv- 
alent interference pattern of 450 lines/mm. Using these values 
and adjusting the derivative of the index of refraction produces 
Fig. 5. 

6. Coordinates by table lookup 

There are numerous instances in which a physical coordinate 
is well defined at each pixel along an image axis, but the re- 
lationship of the coordinate values between pixels cannot be 
described by a simple functional form. An obvious example 
of this would be a three-dimensional image consisting of a se- 
quence of two-dimensional images of an astronomical object 
recorded at an arbitrary sequence of times determined in part by 
weather and time assignment committees. As another example, 
the calibration of some spectrographs, such as those employing 
a diode array detector, is represented best by a list of frequen- 
cies or wavelengths for each pixel on the spectral axis rather 
than some functional form. 

6.1. xxxx-TAB non-linear algorithm 

To support such representations for primary images and im- 
age extensions, we define a table-lookup form for the value of 



CTYPE/a as xxxx-TAB, where xxxx are four letters representing 
the type of coordinate, e.g. TIME or FREQ. As in Paper I, which 
established the "4-3" convention for CTYPE;, the coordinate 
xox is a "real" coordinate, such as FREQ, not an intermediate 
coordinate, such as FREQ-F2W, requiring an additional linear 
or non-linear algorithm in order to be converted into a physical 
coordinate. 



6.1 .1 . -TAB indexing concepts 

Consider first the case of a single (one-dimensional) coordi- 
nate axis. The -TAB algorithm uses a list of coordinate values, 
the coordinate array, to record coordinate values of the appro- 
priate type for the coordinate axis. A second list of the same 
length, the indexing vector, may be used in addressing the co- 
ordinate array. The indexing vector provides one level of indi- 
rection which may be used to vary the sampling frequency of 
the coordinate array along the coordinate axis. The coordinate 
could then be sampled at smaller intervals over that portion of 
the range in which the instrument is more non-linear and sam- 
pled more coarsely over regions in which it is better behaved. 
The indexing algorithm, based on linear interpolation, is de- 
fined in Sect. 6.1.2. If the indexing vector is absent, it is taken 
to have values 1,2,3, ... ,K, where K is the number of values 
in the coordinate array. See Sect. 6.2.2 for an example of the 
use of the indexing vector. 

The concept described above covers separable (one- 
dimensional) axes only. It may be extended simply to M non- 
separable axes so long as the indexing vectors for each of the M 
axes are separable. M coordinate values are required for each 
of the possible index positions. Therefore, the coordinates will 
be in a single (1 + M)-dimensional array. This coordinate ar- 
ray will have dimensions (M, K\,K2, ■ ■ ■ K M ), where K m is the 
maximum value of the index on axis m + 1 of the coordinate 
array. For simplicity, degenerate axes are forbidden; therefore, 
K m > 1. The indexing vectors for each of the M axes each 
contain a one-dimensional array of length K m . 

The data in the indexing vectors must be monotonically in- 
creasing or decreasing, although two adjacent index values in 
the vector may have the same value. However, it is not valid 
for an index value to appear more than twice in an index vec- 
tor, nor for an index value to be repeated at the start or at the 
end of an index vector. Furthermore, repeated index values are 
allowed only in the case of one-dimensional separable axes. 
(See the following section for a discussion of the reasons for 
these limitations and how interpolation is done in such cases; 
see also Sect. 6.2.3 for an example of equal index values and 
further discussion of interpolation.) Application programs must 
not sort the index and coordinate arrays since this makes the 
relative order of the two equal index values indeterminate. The 
requirement for monotonic index values should eliminate any 
need for sorting. Note that it does not imply any ordering of the 
values in the coordinate array. 
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-TAB without an indexing vector 



(fr m = CRVALi + 2 CDLj(pj - CRPIX 7) 



single-row table PSLO 
I 

ip m is a direct index into table cell coordinate array 

jT -1 PS LI 

Ci C2 C3 C4 C5 Cf, C7 ... Cs: 



-TAB with an indexing vector 



i// m = CRVAU + CDELTi £ W-jiPj ~ CRPIX j) 



single-row table PSLO 
I 

Find index t m by interpolating !/(„, in table cell indexing vector "Fj 

PSL2 j, 1 PS LI 

\V l V 1 *Vi'¥^7v ( , y ¥i ■■■ ^ll C i C 2 C 3C 4 C 5 C 6 C 7 ... C K \ 

T,„ selects point in table j 

cell coordinate array Ct„, 



Fig. 6. -TAB logic flow with and without an indexing vector. The coor- 
dinate array subscript m associated with intermediate world coordinate 
axis ( is specified with keyword PVL3. In the case of an independent 
-TAB axis it would have value 1. 

Note that ip m is computed either with the CD Lj form of the linear trans- 
formation matrix or the PC Lj plus CDELT form. 

6.1 .2. Computing -TAB coordinate values 

The indexing algorithm is illustrated schematically in Fig. 6 
for the one-dimensional case. In the general case, to determine 
the M non-separable coordinate values C m , one first determines 
the M indices, ij/ m for addressing the appropriate location in the 
table. If intermediate world coordinate axis i is associated with 
the m th axis in the coordinate array, one evaluates 



i// m = Xi + CRVAL ia, 



(87) 



where x, is computed following the prescriptions of Eq. (1). 
Using linear interpolation, if necessary, in the indexing vector 
for intermediate world coordinate axis i, one determines the 
location, Y m , corresponding to if/ m . Then the coordinate value, 
C m , of type specified by the first four characters of CTYPE/a, 
is that at location (m, Ti, T2, . . . T M ) in the coordinate array, 
again using linear interpolation as needed. 

In detail, the algorithm for computing Y m , and thence C m , is 
as follows. Scan the indexing vector, Q¥ it ¥2, . . .), sequentially 
starting from the first element, *Pi, until a successive pair of 
index values is found that encompass ifr m (i.e. such that < 
ty m < ^k+i for monotonically increasing index values or > 

> *Pj:+i for monotonically decreasing index values for some 
k). Then, when + "P/t+i, interpolate linearly on the indices 
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Vk + l - Vk 



(88) 



which yields an index into the coordinate array. However, if 
^k = y Vk+i(= tym) then the result is undefined. 

In the case where an index value is equal to ifr m , the algo- 
rithm above will find the interval < i// m = *P/t+i for mono- 



tonically increasing index values or >\p m - ^k+i for mono- 
tonically decreasing index values, except when tf/ m = *Pi. Since 
two consecutive index values may be equal, the index follow- 
ing the matched index must be examined. If ^k+i = ^k+i = tym 
(or ¥2 = *Pi = tym), then T m and C m are undefined. 

Linear interpolation via Eq. (88) applies for if/ m in the range 
^1 to inclusive. Outside this range, for K > 1, linear extrap- 
olation is allowed for values of \p m such that v Fi-( v F2- v Pi)/2 < 
\f/ m < »Pi or < \f/ m < X ¥ K + C^k - ^-i)/2 for monotonic 
increasing index values, and for + (¥1 - v F2)/2 > <p m > 
or > if/ m ^ - C^k-i - v I'a:)/2 for monotonic decreasing 
index values. Extrapolation is also allowed for K = 1 with ip m 
in the range Ti - 0.5 < \p m < *¥\ + 0.5 (noting that *Fi should 
be equal to 1 in this case) whence T m = ifr m . 

The value of T m derived from if/ m must lie in the range 0.5 < 
T m < K + 0.5. These extrapolation limits permit assignment of 
coordinates to the regions of the pixels on the boundary of the 
array which are outside of the centers of the boundary pixels 
but within the conceptual "edge" of the boundary pixels. In the 
case of a single separable coordinate with 1 < k < Y„, < k+ 1 < 
K, the coordinate value is given by 



C m - Ck + (T m - k) {Ck+i - Ck) ■ 



(89) 



For T m such that 0.5 < T m < 1 or K < T m < K + 0.5 linear 
extrapolation is permitted, with C m = C\ when K = 1. 

Conceptually, to compute the change in coordinate value 
between and T^+i, in the case when ^+2 = ^k+i and/or 
= Ti, difference the two values of C m obtained for 
<Am = *P/t+i - £ and ty m = ^ + e in the limit that e goes to 
zero. In practice, the computation is straightforward and it is 
not necessary to take limits. 

The inverse computation, in which one determines a \ff m 
given a coordinate C,„, is relatively straightforward, at least in 
the case of a single separable coordinate. One scans the coor- 
dinate vector from the start looking for a pair of values that 
encompass C m . Having found a pair, one must then check the 
indexing vector. If the two index values are unequal, then the 
correct pair has been found. Otherwise the search must be con- 
tinued. When a correct coordinate pair has been found, the in- 
verse of Eq. (89) is applied to determine Y m . The inverse of 
Eq. (88) then yields if/ m . 

It is understood that the general keywords CUNITma, 
CRDERma, and CSYERma apply to the output coordinate C m 
rather than the -TAB coordinate keywords CRVAL ia et al. 
Similarly, if the table lookup determines celestial coordinates, 
the general keywords RADESYSa and EQUINOXa apply to the 
output celestial coordinates rather than the input -TAB coordi- 
nates. 

6.1.3. -TAB implementation, parameters, and 
requirements 

Standard FITS binary tables extensions (XTENSION = 
'BINTABLE') will be used to hold the coordinate array in 
a single cell of a column of a one-row table. The length of 
this array is given in the FITS table header by the repeat 
count in the TFORMn keyword, where n is the number of the 
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Table 7. Parameter keywords used for the -TAB algorithm. 



Keyword 


Use 


Default 


PS«_8a 


table extension name (EXTNAME) 




PVLla 


table version number (EXTVER) 


1 


PVi'_2a 


table level number (EXTLEVEL) 


1 


PSLla 


column name for the coordinate 






array (TTYPEni) 


_ 


PSi'_2a 


column name for the indexing 






vector (TTYPErc 2 ) 


no index 


PVL3a 


axis number (m) in array PSLla 


1 



column containing the coordinate array. The dimensions of the 
coordinate array will be given in the FITS table header by the 
keyword TDIMn set to ' (M , K\ , K2 , ... K M ~) ' , where n is 
the column containing the coordinate array. Note in particular 
that if M = 1 the coordinate array is formally two-dimensional 
though the first axis is degenerate. The repeat count in the 
TFORMn keyword value is the product of M and all the K m . 
The indexing vectors for each of the M axes, if present, will 
occupy separate columns, each containing a one-dimensional 
array of length K m . 

The BINTABLE extension containing the coordinate array 
must be in the same FITS file as the data that reference it. 

The parameters required by -TAB are the table extension 
name (EXTNAME), the table version number (EXTVER), the ta- 
ble level number (EXTLEVEL), the column name for the co- 
ordinate array (TTYPEni), the column name for the indexing 
vector (TTYPE«2), and the axis number m associated with in- 
termediate world coordinate axis i in the coordinate array. The 
keywords used for this purpose are PSLQa, PVLla, PVL2a, 
PSLla, PSL2a, and PVL3a, respectively. These are summa- 
rized in Table 7. For images, PS LQa has no default; for tables a 
missing or blank PS LQa is taken to be the current table (see be- 
low). PVLla and PVL2a both have a default value of 1. PSLla 
never has a default; it must be present and be assigned a value 
actually occurring in table PSLQa. If PSL2a is missing or has 
a value of all blanks, the indexing vector is taken to be a list of 
integers from 1 to K m and if/ m becomes a direct index of axis 
PVL3a in the array specified by PSLla. If PVL2a is present 
and assigned a value, then that value must actually occur in ta- 
ble PS LOa. Note that the values given to PS iAa and PS L2a are 
not case sensitive since the FITS Standard (Hanisch et al. 2001) 
states that "String comparisons with the values of TTYPEn key- 
words should not be case sensitive." 

The use of -TAB for M related axes requires the header 
to specify the array PSLla to be the same for each of the M 
axes. The dimensions of this array must be given in the header 
as (M, K t ,K2, . . . , K M )- These dimensions determine the max- 
imum range of the array index (1 to K m ) for each axis. The 
indexing vectors for each of the M axes, if present, will occupy 
separate columns, each containing a one-dimensional array of 
length K m . The M values of PVL3a must account for all M 
axes. If any of these conditions are not met, the result is unde- 
fined. 



If a FITS file contains multiple XTENSION HDUs (header- 
data units) with the specified EXTNAME, EXTLEVEL, and 
EXTVER, then the result of the WCS table lookup is undefined. 
If the specified FITS BINTABLE contains no column, or multi- 
ple columns, with the specified TTYPEn, then the result of the 
WCS table lookup is undefined. The specified FITS BINTABLE 
must contain only one row. 

No units conversions are to be performed. CUNIT/ must be 
the same as TUNITn of the binary table, where n is the column 
number corresponding to PSLla. 

We recommend strongly that the value chosen for EXTNAME 
always begin with the four letters WCS-. If this is done, generic 
programs will recognize the table as part of the WCS portion 
of the data and will be less likely to separate the table from the 
rest of the WCS. 

The -TAB implementation is more complicated than most 
other WCS conventions because the coordinate system is not 
completely defined by keywords in a single FITS header. 
Software that supports -TAB must be able to gather all the nec- 
essary WCS parameters that are in general distributed over two 
FITS HDUs and in the body of the WCS extension table. The 
producers of FITS data products should consider the capabili- 
ties of the likely recipients of their files when deciding whether 
or not to use the -TAB convention, and in general should use 
it only in cases where other simpler WCS conventions are not 
adequate. 

6.1 .4. -TAB usage in tables 

Binary table extensions containing array data columns may 
need a table-lookup function for coordinate values. Seemingly 
the most convenient form of -TAB would be one in which the 
-TAB array(s) are also table cells in the same row as the data 
array. However, a separate coordinate table would be more eco- 
nomical if it applied to data arrays in multiple rows. The -TAB 
keywords for such tables are /CTYPwa, /Sn_0a, iVnAa, iVnJZa, 
iSnAa, and ;'Sn_2a, where n is the column number of the array 
of data and ; is the intermediate world coordinate axis number. 
If /Sn_0a is missing or blank, the present binary table is taken 
to be the coordinate lookup table. In that case, the coordinate 
array(s) are taken to be single-cell arrays in the same row as the 
data array and keywords iVnAa and NnJZa are ignored. 

Strictly speaking, the -TAB representation is not required 
for binary or ASCII table extensions containing only one data 
value per cell because each of the coordinate values associated 
with the datum may be stored in separate columns. However, 
-TAB does provide a convenient method of solving one of the 
problems pertaining to such tables, namely identifying the data 
column. If the table contains a column of data and another 
one or more columns of coordinate values pertaining to those 
data, then one can define the coordinates of the data column 
as being of type -TAB. In this case, the critical keywords are 
iCTYPna to declare the coordinate axis type and ;S«_la to iden- 
tify the corresponding coordinate column. Since the data value 
column contains only one value, the coordinate column should 
only contain one value. Then, the usual coordinate keywords 
(j'CRPXn / jCRPna, ijPCna, ijCDna, i'CDLTn / iCDEna) may be 
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omitted since their defaults yield the desired result that out- 
put index pixel equals the input data pixel (both of which are 1 
in the assumed case). Note too that tables of this type — one 
value per table cell — may be in either ASCII or binary table 
form. 



6.2. -TAB examples 

To illustrate the use of -TAB with two non-separable axes, let us 
consider a specific example of a four-dimensional image array. 
Assume that axes 1 and 3 are handled either linearly or by one 
of the non-linear single-axis cases (including -TAB). However, 
assume that axes 2 and 4 require table lookup in a mutually 
dependent fashion. Thus ;' = 2 for m = 1 and i = 4 for m = 2 
and the required keywords are 
PS2JJ = 'WCS-TAB' PS4JS = 'WCS-TAB' 

PV2_1 = 1 PV4_1 = 1 

PV2_2 = 1 PV4_2 = 1 

PS2_1 = 'COORDS' PS4_1 = 'COORDS' 

PS2_2 = ' INDEX2 ' PS4_2 = ' INDEX4 ' 

PV2_3 = 1 PV4_3 = 2 

where the first four keyword pairs must match exactly and the 
last two pairs must have different values. 

A real example that could use this particular algorithm is 
represented by a spectral image of a portion of the sky taken 
with a single radio telescope. The observer commands the tele- 
scope to point at a regular grid of coordinates, but, due to wind 
loading and other pointing errors, the telescope achieves the 
commanded positions only approximately. The actual celestial 
coordinates observed are, however, accurately measured. These 
data could then be represented usefully as a two-dimensional 
array of spectra, but accurate celestial coordinates for each 
spectrum could be found only by a 2-dimensional table lookup. 

Since there has been no generally agreed upon FITS for- 
mat for spectral data with explicit wavelengths assigned to each 
pixel, data providers have resorted to defining their own for- 
mats. Examples from the Hubble Space Telescope and the Very 
Large Array are shown below, recast into the -TAB algorithm. 

6.2.1. -TAB examples: HSTdata 

Two types of Hubble Space Telescope data serve as examples 
of the -TAB algorithm. The two cases illustrate the evolution 
from simple images to the more powerful constructs provided 
by FITS extensions and binary tables. The purpose of dis- 
cussing these formats is not to explain the HST formats, but 
to illustrate a couple of types of data that can be represented by 
the -TAB algorithm. Therefore, some details of these formats 
are ignored. 

The early HST spectrographs, FOS and GHRS, adopted a 
format based only on simple FITS images. The basic concept 
is that the spectral flux values are given in one image and the 
vacuum wavelengths, in Angstroms, are given as the pixel val- 
ues in another image. The two are associated by filenames. The 
filenames have the same basename, but different filename ex- 
tensions for an exposure. To represent these spectra with a FITS 
WCS based on the -TAB method, while continuing to use an 



Table 8. WCS keywords in the spectral image for FOS and GHRS 
example. 

CTYPE1 = 'WAVE-TAB' / Coordinate type 

CUNIT1 = 'Angstrom' / Coordinate units 

PS1J8 = 'WCS-TAB' / Coord table extension name 

PS1_1 = 'WAVELENGTH' / Coord table column name 

Table 9. WCS keywords in the STIS table for two of the spectral 
columns. 

1CTYP4 = 'WAVE-TAB' / Coordinate type 

1CUNI4 = 'Angstrom' / Coordinate units 

1S4_1 = ' WAVELENGTH ' / Coords for gross spectrum 

1CTYP5 = 'WAVE-TAB' / Coordinate type 

1CUNI5 = 'Angstrom' / Coordinate units 

1S5_1 = 'WAVELENGTH' / Coords background spectrum 



image representation for the spectral fluxes, one replaces the 
separate wavelength image with a table extension. 

Table 8 shows the minimum WCS keywords required in 
the primary image header. The PS1_8 keyword value defines 
the coordinate table to be in a binary table extension with name 
' WCS-TAB ' . The coordinate table extension consists of one col- 
umn, named 'WAVELENGTH', and one row. The array values 
and length correspond to the original image format. 

Other WCS keywords defined in Paper I would also be in- 
cluded as needed. The intermediate coordinate transformation 
must produce indexing values equal to the pixel coordinate in 
the image. Since the defaults for these keywords are defined in 
Paper I to produce pixel coordinates, it is not strictly necessary 
to include these keywords. 

The second generation HST spectrograph, STIS, adopted a 
binary table format very close to one of those defined for the 
-TAB algorithm. In the STIS one-dimensional extracted echelle 
format each exposure is stored in a separate binary table. The 
extracted data from an exposure consist of a number of echelle 
orders each of which is stored in a separate row. The wave- 
lengths and various associated "spectra" are array elements in 
different columns. An associated spectrum is an array of a sin- 
gle type data quantity such as fluxes, errors, or data quality 
flags. 

To convert this format to the -TAB representation only re- 
quires adding the ICTYPn and lSn_l keywords to specify the 
coordinate type and coordinate column name for each spectral 
column n. There are six types of spectral quantities for each 
order: gross, background, net, flux, error, and data quality flag. 
So there must be six pairs of keywords. Since the same wave- 
length column applies to each of the spectra, the keyword val- 
ues are repeated six times, but with different column numbers 
in the keywords. Table 9 shows the WCS keywords for two of 
the spectral columns. The 1S«_0 keywords are omitted, signi- 
fying that the coordinate array column is in the same table as 
the spectra. 

As noted for the FOS and GHRS format, additional WCS 
keywords are included as needed and the intermediate coordi- 
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Table 10. Sample coordinate table. Each displayed column is actually 
a one-cell array of the (single-row) table. 
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nate transformation values may be absent since this defines an 
indexing by pixel coordinate. In the STIS data, the units of the 
wavelength array are vacuum Angstroms so the lCUNIn key- 
word would be required. 

6.2.2. -TAB example: radio interferometry 

Radio interferometry data are now best represented in binary 
tables with each row representing the visibility at one time on 
one antenna pair with a vector of data representing the complex 
values for all polarizations and frequencies observed. These 
telescopes often allow the user to observe N( regularly-sampled 
spectral channels starting from L arbitrarily chosen interme- 
diate frequencies (IFs). Furthermore, this spectral data pattern 
may be used with one receiver and, a few minutes later, with 
other receivers at very different frequencies. Because of the 
great flexibility in the choice of the IFs and receivers, the fre- 
quencies present in the data may only be described by a table. 
Because that description is very repetitive, we choose to put it 
in a table separate from the visibility table. 

We can construct a single column table of all Nt fre- 
quencies observed. However, by using a second column for an 
indexing vector, we can take advantage of the regular sam- 
pling of the N( spectral channels. If we set i'CRVLn, i'CDLTn, 
and /CRPXn (for binary table format) or CRVAL;, CDELT/, and 
CRPIX/ (for random groups and image formats) to one, then 
Table 10 may represent the actual frequencies, where V( is the 
frequency of channel 1 in IF I and 5c is the increment in fre- 
quency between the regularly sampled channels of IF I. This 
example is shown graphically in Fig. 7. 

6.2.3. -TAB examples: multiple exposures and 
sampling 

In general, the FITS WCS papers do not consider the pro- 
cess by which the physical effects represented by the values 
in the FITS array have been measured. The measurable quan- 
tity may have been sampled over regions small enough to be 
considered as points in a widely separated grid of world coor- 
dinates. Alternatively, the measurable quantity may have been 
integrated across overlapping adjacent regions. The FITS WCS 
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-TAB parameters and values 

CRVALi = CRPIX! = CDELT; = PC(_( = 1 

single-row table PSLO 
values in the cell of the indexing vector PS i_2 = 

|1, 7, 8, 11, 12, 18, 19, 25, 26, 30 1 

corresponding values 
in the cell of the coordinate array PS LI = 

vi, vi + 65i, v 2 , v 2 + 3<5 2 , v 3 , v 3 + 65 3 , v 4 , v 4 + 6S4, v 5 , v 5 + 4i5 5 



Fig. 7. Example taken from radio interferometry using -TAB with an 
indexing vector. The FITS keywords shown are suitable for the ran- 
dom groups format. The observation is made at a number of frequen- 
cies, with a number of regularly spaced St spectral channels begin- 
ning from each of a number of arbitrary base frequencies v e . The use 
of an indexing vector reduces the number of values in the table from 
one array of 30 to two arrays of 10 each. In a real case the number of 
spectral channels would be significantly larger, making the space sav- 
ings significant. In this example, pixel pj = 6 produces \j/ m = 6. This 
lies at T„, = l|. The resulting coordinate is then ^vi + | (vi + 6<5i) 
which, as one would expect, equals v ( + 5<5i . Note that this example 
involves only a single independent -TAB axis, so that PVL3 = m = 1. 



papers do not indicate whether it is valid to presume that the 
values in adjacent array elements represent physically adjacent 
quantities, nor do they indicate whether it is valid to interpolate 
the WCS across the extent of an array element. 

Some world coordinate axes are defined to have integral 
values with conventional meanings (e.g. COMPLEX, STOKES). 
Interpolating a world coordinate value across the extent of a 
FITS array element along the direction of an axis which is con- 
ventionally defined to be integral would clearly be inappropri- 
ate. Nevertheless, for other types of world coordinate axes in 
traditional FITS arrays, the adjacent elements in the NAXIS- 
dimensional array can often be presumed to represent adjacent 
locations in a measurement continuum. This adjacency is rele- 
vant, for example, to image display programs when the image is 
magnified so that one FITS array element extends across mul- 
tiple display pixels. In such cases a FITS display program may 
attempt to interpolate WCS values across the extent of an ar- 
ray element. With CTYPEA: values other than -TAB this is com- 
monly done by linear interpolation between the presumably 
adjacent array elements. However for coordinate axes where 
the CTYPEfe is described by the xcxc-TAB non-linear algorithm 
there can be no presumption of adjacency. 

Notwithstanding these caveats, the use of -TAB does per- 
mit specification of changes in the coordinate value across the 
extent of a single array element. Consider a sequence of two- 
dimensional images of the sky which have been re-gridded 
such that the pixel array for each image shares the same 2-D 
WCS for celestial coordinates. A single 3-D FITS array can 
hold such a sequence of 2-D images. If the same instrument 
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Table 11. Top: keywords in main FITS header for -TAB example, bottom: table keywords and sample arrays for that example. 
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v. 1 lrLj 




'WAVF-TAR' 
minvL 1 /ID 




/ 


Qnpp+'p^l avi c l*iv ■fsVilp 1 not— nn 
o^Jtr l. lx d j. a. a x z> uy LcLUJ.tr iuujv u.jj 


CTYPE4 




' TTMF-TAB ' 




/ 
/ 


1" Pitir>fir3 1 ayi c l*iv taVilp 1 not— nn 
L triLiuui di ctA J. s uy LcLUJ.tr iuujv u.|J 






' RanHna ' 

LJ CL1 L L I I J CL kD Z> 








CNAME4 




' Observation Date ' 








CUNIT3 




' m ' 




/ 


Wavcl ferity Lilo XI L ILLtr L trX s 


CUNIT4 








/ 


oVi con?a 1" i on H 3 f" p c "in vp a T 1 c 
UUr3trJ.Va.LJ.UIl UdLtrr> X1L ytrdXra 


CRPIX3 







_ 5 


/ 


Qnpp+'p^il hin itii ti al pHnp n~F initial niirpl 

JS^Jtr L. LX dX UJ.1L XlLXLXdX tr U.LJ tr Ul J-lLJ.LJ.ClJ- X AC X 


CRPIX4 






i 
x 


/ 


LtritipUI cL-L UXILb di tr a. L UtrLJJ.IL/ trlLU. LXHltr 


CDELT3 






i 
i 


/ 
/ 


r>ptrCHI dX X 1 LLltr A dl X dy X r> Ulltr lull L ptrX UXll 


CDELT4 






I 


/ 


+• rn nm" 3 1 "inHpY arrav i c nn p i in i +■ riPT 1 hi n 
LtrlLLUUX dX X 1 LU.tr A dX X dy X Ul Ltr UX LXL JJtr X UX1L 


CRVAL3 







. j 


/ 


r>ptrL. LX dX X trltrX trILL.tr Is L.trILLtrl UI X dUXU trApUraUX tr 


CRVAL4 






fi) 
\y 


/ 


tDTTiiinral rflfflraTiro i c ctart r»"F raHi n PYnriCimp 
LtrlLLUUX dX XtrXtrXtrlLL.tr Xr5 b Lai L UX X dUXU trAUUraUXtr 








1 
1 


/ 
/ 


jia pixel axis increments sra image coord axis 


pr^ 4 






8 






pr4 3 






1 


/ 


3rd pixel axis increments 4th image coord axis 


pr4 4 






1 


/ 


4th pixel axis is degenerate (1 point) 


PQ 3 Ifl 




IfllLi-LdDie 








PS3_1 




'WaveCoord' 








PS3_2 




'Wavelndex' 








PS4JS 




'WCS-table' 








PS4_1 




'TimeCoord' 








PS4_2 




'Timelndex' 









EXTNAME = 'WCS-table' 
EXTVER = 
EXTLEVEL= 



TF0RM1 = '8E' 
TTYPE1 = 'Wavelndex' 



(8-5, 
1-5, 
1-5, 
2.5, 
2.5, 
3.5, 
3.5, 
4.5) 



1 
1 

TUNIT2 = 'm' 
TF0RM2 = '8D' 
TTYPE2 = 'WaveCoord' 



(0.21106114, 
0.21076437, 
2.8E-6, 
2.2E-6, 
500.0E-9, 
650.0E-9, 
1.24E-9, 
2.48E-9) 



TF0RM3 = '8E' 
TTYPE3 = 'Timelndex' 



(0.0, 
1.0, 
1.0, 
2.0, 
2.0, 
3.0, 
3.0, 
4.0) 



TUNIT4 = 'a' 
TF0RM4 = '8D' 
TTYPE4 = 'TimeCoord' 



(1997.84512, 
1997.84631, 
1993.28451, 
1993.28456, 
2001.59234, 
2001.59239, 
2002.18265, 
2002. 18301) 



produced each of the images on different dates then the third 
dimension might be purely temporal. Alternatively, if differ- 
ent instruments (radio, IR, optical, X-ray) produced each of the 
images then the third dimension might be spectral. In practice, 
however, images in such a spectral 3-D image will typically 
also have different observation dates, so the image will actu- 
ally span four dimensions of world coordinates. 

Akin to the example of a long slit spectrum in Paper II, 
the WCS of this four-dimensional case can be represented by 
the 3-D array. This fourth dimension could be represented 
with NAXIS = 4 and NAXIS4 = 1, or with NAXIS = 3 and 
WCSAXES = 4. Whether the sequence of images is temporal 



or spectral, it is certain that the exposure duration or the spec- 
tral bandpass have a non-zero extent and the -TAB representa- 
tion does provide the ability to communicate the start and end 
times of an observation, or the minimum and maximum wave- 
lengths of an observation. For this example the WCS keywords 
are shown in Table 11, where the keywords relating to axes 1 
and 2 are omitted because they are simply celestial coordinates 
as described by Paper II. In this example, note the clever use 
of the PC matrix to cause the fourth coordinate axis to depend 
on the third pixel axis, but not the degenerate fourth pixel axis. 
FITS writers can use details such as these to communicate how 



1 IJ 
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Table 12. Recognized values for SPECSYSa, SSYSOBSa, and SSYSSRCa. 



Value 


Definition 


Magnitude 


I 


b 


Reference 


TOPOCENT 


Topocentric 


0.0 km s" 1 


- 


- 




GEOCENTR 


Geocentric 


0.5 


- 


- 




BARYCENT 


Barycentric 


30 


- 


- 


Stumpff (1980) 


HELIOCEN 


Heliocentric 


30 


- 


- 


Stumpff (1980) 


LSRK 


Local standard of rest (kinematic) 


20 


t56 


23 


Delhaye (1965), Gordon (1976) 


LSRD 


Local standard of rest (dynamic) 


16.6 


53 


25 


Delhaye (1965) 


GALACTOC 


Galactocentric 


220 


90 





Kerr & Lynden-Bell (1986) 


LOCALGRP 


Local Group 


300 


90 





de Vaucouleurs (1976) 


CMBDIPOL 


Cosmic microwave background dipole 


368 


263.85 


48.25 


Bennett et al. (2003) 


SOURCE 


Source rest frame 


any 









t LSRK is normally quoted as right ascension 18 hours, declination 30 degrees (1900). 



a display program might provide meaningful world coordinates 
at locations across an array element. 4 

In order to provide coordinate values across the entire array 
element, the coordinate values on the boundaries between the 
array elements are multiply defined. Such a result should be 
expected for array elements representing non-adjacent physi- 
cal values. It is evident in this case that the order of the index 
arrays should not be permuted during the coordinate lookup. 
The lower portion of Table 1 1 shows a subset of the content of 
the table(s) referred to by PS 3_0 and PS4_6>. The four columns 
in the figure are meant to illustrate the association of table 
header keywords and table values. In the FITS file, the key- 
words would appear in the table header along with many other 
keywords while the four arrays of values would appear in four 
cells of the one-row binary table. 

The indexing vectors in the table of Table 1 1 illustrate the 
reason to allow two values within the indexing vector to be 
identical. The rules for the linear interpolation within such vec- 
tors are clear. One must select the two index locations having 
values immediately surrounding i// m . Thus in the present exam- 
ple, if the time coordinate has tf/ m = 1.1, the third and fourth 
index and coordinate values are interpolated. The output time 
coordinate is 

1993.28451 + ^"n I i'o ( 1993 - 28456 - 1993.28451)- 

The result of the table look up is undefined if the ifi m is equal to 
the repeated value in the indexing vector. 

7. Coordinate reference frames 

Frequencies, wavelengths, and apparent radial velocities are 
always referred to some selected standard of rest (reference 
frame), and while they are measured, of necessity, in the ob- 
server's rest frame, they may also be corrected to some other 
frame. 5 The velocity correction is computed from the vector 

4 The precise interpretation of CTYPEfc keywords that describe TIME 
is deferred to a possible future paper on temporal coordinates in FITS. 

5 Proper relativistic formulas should be used; see Gibbs & Tao 
(1997) or Rindler (1982). Better yet, observers should consider using 



dot product of the direction vector and the relative velocity 
vector of the two reference frames. In addition, the corrections 
from topocentric, the frame in which the observations are usu- 
ally made, to geocentric and then to barycentric or heliocentric 
are dependent on the dot product with time-variable velocity 
vectors. As a consequence, the "velocity" with respect to the 
reference frame depends on direction. Differential effects over 
a field of view may be important in some high precision ap- 
plications; for example, they may amount to 5 km s _1 for the 
Local Group correction over a field of one degree. In another 
example, observations of Galactic CO over two-degree fields 
separated by two degrees failed to align by two spectral chan- 
nels (Mangum et al. 2000). 

If a frequency, wavelength or apparent radial velocity asso- 
ciated with an image plane has been corrected to some other 
frame by transforming the CRVALfca et al. values, then differ- 
ential errors may arise at points away from the reference point. 
For example, it is normal in radio astronomy to observe a field 
while holding constant the velocity of the reference point with 
respect to the local standard of rest. See Sect. 10.2 for a more 
detailed discussion of the complications inherent in such data. 

Nonetheless, each image plane shares a constant topocen- 
tric frequency (or apparent radial velocity). The velocity {and 
frequency and wavelength) with respect to the local standard 
of rest is then a function of celestial coordinate within the im- 
age. In order to denote this, we introduce two character-valued 
keywords. The first, 

SPECSYSa (character-valued) 

describes the reference frame in use for the spectral-axis coor- 
dinate^). The second, 

SSYSOBSa (character-valued) 

describes the spectral reference frame that is constant over 
the range of the non-spectral world coordinates and has de- 
fault value ' TOPOCENT ' . The recognized values are given in 

standard software described by Lindegren and Dravins (2003) to per- 
form the conversion from topocentric coordinate to the standardized 
Barycentric Celestial Reference System. 
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Table 12. In this table, the magnitude column gives the approx- 
imate magnitude of the velocity vector that defines the particu- 
lar reference frame and the I and b columns give the standard 
Galactic longitude and latitude of the reference frame, Blaauw 
et al (1960). This table includes the traditional heliocentric sys- 
tem which the authors believe should be deprecated. It has fre- 
quently been misused as an alias for barycentric and would, if 
high accuracy is desired, require the observing date in order 
to allow conversions between it and the other systems. The ta- 
ble also includes the WMAP microwave background dipole as 
a velocity system. Since, as the data are refined over time, this 
system will change, its use at present should be restricted to sit- 
uations for which it is particularly appropriate and the parame- 
ters used should be well commented. The radio-astronomy ex- 
ample described above has SPECSYS = 'LSRK' and SSYSOBS 
= ' TOPOCENT ' indicating that the spectral axis is expressed as 
a kinematic Local Standard of Rest velocity, but that each right 
ascension by declination image plane is at a constant topocen- 
tric spectral coordinate. 

Many spectrometers suffer from a significant variation in 
their spectral coordinate as a function of celestial coordinate. 
Such variations are the primary subject of Paper IV, while the 
concepts of this work refer to images produced by, or corrected 
to, an ideal spectrometer. In general we recommend that the 
frame corrections will get folded into the instrumental ones 
and that the image produced after the Paper IV corrections will 
have identical SSYSOBS and SPECSYS. This is not required, but 
would certainly simplify matters. 

It is not appropriate for this work to define the parameters 
of the recognized rest frames. The parameters needed to com- 
pute geocentric frequencies/velocities from topocentric are the 
sidereal time (or Earth rotation angle) and the observatory lo- 
cation. The observing date is needed to convert from geocentric 
to barycentric coordinates. The new keywords 

MJD-AVG (floating-valued) 

and 

DATE-AVG (character- valued) 

are reserved to give a representative time for the whole obser- 
vation, suitable for calculating proper motion, apparent place, 
local apparent sidereal time, velocity with respect to the center 
of the Earth and barycenter, etc. The DATE-AVG keyword shall 
be expressed in the manner described by Hanisch et al. (2001) 
for the DATE-OBS keyword. If both keywords are present, they 
must agree to adequate accuracy or the result is undefined. If 
both keywords are absent, DATE-OBS may be used to determine 
the representative time, although DATE-OBS refers to the begin- 
ning of an observation and hence is not entirely suitable. For 
high-precision applications, a time system to apply to MJD-AVG 
must be defined in the manner proposed by Bunclark & Rots 
(1997). Note also that, for long integrations, no single time is 
adequate for computing proper motion, sidereal time, and the 
like. 

Astronomers have traditionally quoted telescope positions 
using terrestrial longitude, latitude, and height, but, for com- 
putation of topocentric velocity, these values are immediately 



converted into a geocentric Cartesian triple. A precise con- 
version from longitude, latitude, and height, however, requires 
knowledge of the geodetic datum in which they are expressed. 
There are about 1000 geodetic datums in use. These are de- 
scribed by the United States Department of Defense (2000) in 
a printed document also available over the Internet. Some of 
the datums deviate from geocentricity by over a kilometer. 

Rather than burden FITS with these details of geodetic tra- 
dition, three new keywords 

OBSGEO-X (floating-valued) 
0BSGE0-Y (floating-valued) 
0BSGE0-Z (floating-valued) 

are reserved to define a representative location for the observa- 
tion in a standard terrestrial reference frame. The location val- 
ues should be right-handed, geocentric, Cartesian coordinates 
valid at the epoch of MJD-AVG or DATE-AVG. Position errors of 
several kilometers should have negligible impact on the com- 
putation of the diurnal velocity correction, so for this purpose 
geodetic accuracy is not a requirement. However, it should be 
possible to provide coordinates with an accuracy of lm or bet- 
ter based on a recent satellite-based geodetic reference frame 
and we recommend that FITS writers do so. Although it is be- 
yond the scope of this paper to define particular terrestrial ref- 
erence frames or tectonic models, only frames based on the 
ITRS (McCarthy & Petit 2004) are suitable for high-precision 
applications. Most post-1980 reference frames, including GPS, 
agree with recent versions of the ITRF to about 0. 1 m. Web 
sites 

http : //earth- info . nga .mil/GandG/ 
http : //www . ng s . noaa . gov/TOOLS/XYZ/xyz . html 
provide general geodetic information and a solution based on 
the outputs of almost any hand-held GPS unit, respectively. 
These references, together with information about the origin 
of the traditional telescope position, should allow the determi- 
nation of the keyword values defined in this paper. 

It may be helpful for the FITS writer to provide information 
on the relative radial velocity between the observer and the se- 
lected standard of rest in the direction of the celestial reference 
coordinate. The keyword 

VELOSYSa (floating-valued) 

is hereby reserved for this purpose; its units shall be ms~'. 
CUNITfea is not used since the WCS version a might not be 
expressed in velocity units. 

The frame of rest defined with respect to the source 
(SOURCE) is useful for comparing internal apparent radial ve- 
locities in objects at different redshifts. This allows the FITS 
writer to apply all the cosmological and other model-dependent 
corrections, leaving a coordinate local to the object. For this 
frame of rest, it is necessary to define the velocity with respect 
to some other frame of rest. The keywords 

ZSOURCEa (floating-valued) 

and 

SSYSSRCa (character- valued) 



Table 13. New spectral keywords including forms for use in tables. All table keywords except the coordinate name have the same form for the 
BINTABLE vector image format and the pixel list format. 



Primary 


BINTABLE Pixel 


Type 


Sect. 


Keyword 


Array 


list 






Description 


fXTAMr if, 
t-lNAllE, I CI 


l^NRllCl 1 {.NAIJCl 


character 




Coordinate axis name 






floating 


1 A A 


Line rest frequency (Hz) 






floating 


1 A A 


Line rest frequency alternate for primary array only 




l\Wri vnu 


floating 


1 A A 


Line rest wavelength in vacuum (m) 


ircLi I it/ 




character 


7 

1 


Spectral reference frame (from Table 12) 


ccvcnDC/. 




character 


-1 

1 


Spectral reference frame fixed during observation (from Table 12) 


rlJU-RM\j 




floating 


/ 


Average date of observation (Julian date— 2400000.5) 


UAL t-Avu 




character 


7 

/ 


Average date/time of observation 


UDibtU-A 




floating 


7 

/ 


Observation X (m) 






floating 


7 


Observation Y (m) 






floating 


7 
/ 


Observation Z (m) 


VELOSYSa 


VSYSrca 


floating 


7 


Radial velocity wrt standard of rest (ms -1 ) 


ZSOURCEa 


ZSOUrca 


floating 


7 


Redshift of source for SOURCE cases (unitless) 


SSYSSRG? 


SSRCrca 


character 


7 


Spectral reference frame for SOURCE cases (from Table 12) 


VELANGLa 


VANGrca 


floating 


A 


Angle of true velocity from tangent to line of sight 



are hereby reserved; they should be used to document the 
adopted value for this systemic velocity of the source. 
ZSOURCEa is a unitless redshift; SSYSSRCa specifies the ref- 
erence frame for this parameter and may have any value 
in Table 12 except SOURCE. The new floating-valued and 
character-valued keywords hereby reserved are listed 6 in 
Table 13 along with the other keywords defined in this paper. 

8. Alternate FITS image representations: pixel list 
and vector column elements 7 

The use of coordinate keywords to describe a multi- 
dimensional vector in a single element of a FITS binary table 
(Cotton et al. 1995), and a tabulated list of pixel coordinates in 
a FITS ASCII (Harten et al. 1988) or binary table was discussed 
in Papers I and II. In this section, we extend those discussions 
to the keywords specific to spectral coordinates. This conven- 
tion is considered an integral part of the full world coordinates 
convention. 

8.1. Keyword naming convention 

Table 13 lists the corresponding set of coordinate system key- 
words for use with each type of FITS image representation. The 
keywords defined in this paper and their allowed values are ap- 

6 The frame and velocity keywords in Table 13 apply to terres- 
trial observers only. Observations from spacecraft require additional 
description of their motions - either a position and velocity or an 
ephemeris. The required keywords and/or tables for this are beyond 
the scope of this paper. 

7 Suggested by contributions from William Pence, Arnold Rots, 
and Lorella Angelini of the NASA Goddard Space Flight Center, 
Greenbelt, MD 20771. 



plicable to all three image storage formats. The following notes 
apply to the naming conventions used in Table 13: 

- a is a 1 -character coordinate version code and may be blank 
(primary) or any single uppercase character from A through 
Z. 

- n is an integer table column number without any leading 
zeros (1-999). 

- ;' is a one-digit integer axis number (1-9). 

When using the BINTABLE vector image format, if the ta- 
ble only contains a single image column or if there are multiple 
image columns but they all have the same value for any of the 
keywords in Table 13 except CNAME/a, then the simpler form 
of the keyword name, as used for primary arrays, may be used. 
For example, if all the images in the table have the same loca- 
tion, then one may use one set of OBSGEO-X, OBSGEO-Y, and 
OBSGEO-Z keywords rather than multiple OBSGXn, OBSGYn, 
and OBSGZn keywords. If both forms of these keywords ap- 
pear, then column-specific values are applied to the specified 
columns and the non-specfic values are applied to all other 
columns. (Note, however, that the WCS keywords defined for 
tables in Papers I and II must always be specified using the 
more complex keyword name with the column number suffix 
and the axis number prefix.) 

9. Spectral coordinate variation with other 
coordinates 

There are instruments having spectral coordinates that are a 
function of other coordinates in the resulting image, e.g. ce- 
lestial position. For example, optical astronomers use imag- 
ing Michelson interferometers and scanning Fabry-Perot in- 
struments to produce three-dimensional images with two celes- 
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tial and one spectroscopic coordinate. Because the path length 
through such instruments increases off the optical axis of the 
telescope, the observed wavelength is a function of celestial 
coordinate in each plane of the image cube. In general, this is a 
subject reserved for Paper IV since such instruments almost al- 
ways have lesser distortions that need to be parameterized and 
corrected along with the path length correction. 

One particularly difficult case should be mentioned here. 
Objective prism instruments produce data in which the image 
of each object in the field has been spread out into a spectrum. 
Thus the output image is the convolution of the sky with both 
a spatial point-spread-function and a spectral dispersive point- 
spread function, both in the two dimensions of the output im- 
age. The methods used to specify coordinates described in the 
present paper and in Papers I and II are unable to handle such 
complex data. However, rather than ignore the world coordi- 
nates entirely, writers of such images should consider provid- 
ing the celestial coordinates that would be correct for the image 
if it had been taken at a single representative wavelength. That 
wavelength would be specified in a third WCS axis. Alternate 
axis descriptions could then be used to specify the same in- 
formation for other wavelengths. Multi-order echelle spectro- 
graphs also produce images in which multiple spectral orders 
overlap each other as well as the celestial axes and, hence, their 
coordinates cannot be described by the methods of Papers II 
and IH. 

Analysis of both objective prism and echelle spectrograph 
images can produce tables of spectra that are easily described 
by the present methods. 

10. Example 

10.1. Basic spectral-line header 

Let us consider an example. The partial FITS header in 
Table 14 was produced for data from the Very Large Array. 
The FITS WCS standard states that all WCS keywords must 
be reproduced for the alternate descriptions, but those for non- 
spectral axes have been omitted for the sake of brevity and 
clarity. We have also updated the keywords to the standards of 
Papers I, II, and III. We will derive additional header keywords 
that could also have been written. These are shown in Table 15, 
again leaving out axes 1 and 2 of the alternate WCS versions. 

The spectral-line cube is regularly sampled in frequency 
and the primary description is of a linear frequency axis (from 
the CTYPE3 keyword value 'FREQ'). CRVAL3 and CRPIX3 
record that the 63-channel spectrum was centered with chan- 
nel 32 at 1.37835117405 GHz, some 42 MHz from the rest 
frequency (RESTFRQ = 1.420405752 GHz) of neutral hydrogen 
(HI), the line under investigation. As indicated by the SPECSYS 
keyword, it is shown as if it had been observed at a constant 
topocentric frequency. For the moment, let us assume that the 
observations were made over a short interval and that this is 
correct. 

As allowed by Paper I, we use alternate axis description 
letter codes in a mnemonic fashion. Using Z for optical, the 
first alternate axis description for the spectral axis reflects the 
observer's request that a Doppler shift be applied to the HI line 



Table 14. Example image header. 



123456789 123456789 123456789 123456789 123456789 



frequency appropriate for a source with a barycentric optical- 
convention velocity of 9120km s -1 . Combining Eqs. (30) and 
(6), we get 

v = v /(l+Z/c), 

whence the barycentric HI rest frequency shifts to 
1.37847122 GHz, some 120 kHz greater than the topocentric 



SIMPLE = T 

BITPIX = -32 

NAXIS = 3 

NAXIS1 = 1024 

NAXIS 2 = 1024 

NAXIS 3 = 63 

EXTEND = T / Tables may follow 

OBJECT = ' 3C353 ' 

TELESC0P= 'VLA' 

DATE-0BS= '1998-09-29' / Start obs 

BUNIT = 'Jy/beam' 

EQUINOX = 2.000000000E+03 

DATAMAX = 5.080558062E-01 

DAT AM IN = -1. 179690361E-01 

WCSNAME = 'TopoFreq' 

CTYPE1 = 'RA— -SIN' 

CRVAL1 = 2.60108333333E+02 

CDELT1 = -2.777777845E-04 

CRPIX1 = 512.0 

CUNIT1 = 'deg' 

CTYPE2 = 'DEC--SIN' 

CRVAL2 = -9.75000000000E-01 

CDELT2 = 2.777777845E-04 

CRPIX2 = 513.0 

CUNIT2 = 'deg' 

CTYPE3 = 'FREQ' 

CRVAL3 = 1.37835117405E+09 

CDELT3 = 9.765625000E+04 

CRPIX3 = 32.0 

CUNIT3 = 'Hz' 

0BSGE0-X= -1.601185365E+06 / [m] 

0BSGE0-Y= -5.041977547E+06 / [m] 

0BSGE0-Z= 3.554875870E+06 / [m] 

MJD-AVG = 51085.979 

SPECSYS = 'TOPOCENT' 

RESTFRQ = 1.428405752E+09 / [Hz] 

RADESYS = ' FK5 ' 



CNAME3Z = 'Barycentric optical velocity' 
RESTWAVZ= 0.211061139 / [m] 

CTYPE3Z = ' V0PT-F2W 
CRVAL3Z = 9.120000E+06 
CDELT3Z = -2. 1882651E+04 

CRPIX3Z = 32.0 

CUNIT3Z = 'm/s' 

SPECSYSZ= 'BARYCENT' 

SSYS0BSZ= 'TOPOCENT' 
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Table 15. Additional keywords for example header. 



123456789 


123456789 123456789 123456789 123456789 


VELOSYSZ= 


26.1O8E+03 / [m/s] 


CNAME3F = 


'Barycentric frequency' 


CTYPE3F = 


' FREQ ' 


CRVAL3F = 


1.37847121643E+89 


CDELT3F = 


9.764775E+04 


CRPIX3F = 


32.0 


CUNIT3F = 


'Hz' 


SPECSYSF= 


' BARYCENT ' 


SSYSOBSF= 


'TOPOCENT' 


VELOSYSF= 


26.108E+03 / [m/s] 


CNAME3W = 


'Barycentric wavelength' 


CTYPE3W = 


'WAVE-F2W 


CRVAL3W = 


9.217481841862 


CDELT3W = 


-1.5405916E-05 


CRPIX3W = 


32.0 


CUNIT3W = 


'm' 


SPECSYSW= 


' BARYCENT ' 


SSYSOBSW= 


'TOPOCENT' 


VELOSYSW= 


26.108E+03 / [m/s] 


RESTFRQR= 


1.420405752E+09 / [Hz] 


CNAME3R = 


'Barycentric radio velocity' 


CTYPE3R = 


' VRAD ' 


CRVAL3R = 


8.85075090419E+06 


CDELT3R = 


-2.0609645E+04 


CRPIX3R = 


32.0 


CUNIT3R = 


'm/s' 


SPECSYSR= 


' BARYCENT ' 


SSYSOBSR= 


'TOPOCENT' 


VELOSYSR= 


26.108E+Q3 / [m/s] 


CNAME3V = 


'Barycentric apparent radial velocity' 


RESTFRQV= 


1.420405752E+09 / [Hz] 


CTYPE3V = 


' VEL0-F2V 


CRVAL3V = 


8.98134229811E+06 


CDELT3V = 


-2.1217551E+04 


CRPIX3V = 


32.0 


CUNIT3V = 


'm/s' 


SPECSYSV= 


'BARYCENT' 


SSYSOBSV= 


'TOPOCENT' 


VELOSYSV= 


26.108E+03 / [m/s] 



frequency given by CRVAL3. The difference, of course, simply 
reflects the instantaneous relative velocity of the barycentric 
and topocentric reference frames at the time of observation. 
Using standard software (the STARLINK program rv), with 
the observatory position defined by OBSGEO-X, Y, Z, at the 
time indicated by MJD-AVG, and for the source at (J2000) 
right ascension and declination 17 h 20 m 26'\ -00°58'30" 
indicated by CRVAL1 and CRVAL2, it may be verified that the 
topocentric correction amounts to 26. 108 km s~'. Now this 
correction between reference frames is a true velocity and thus 
the relativistic formula, Eq. (8), should be used, whereupon 



the frequency recorded by CRVAL3 may be obtained. The 
frequency decreases as the observer moves away from the 
source position. 

With CRVALZ set to 9.12 x 10 6 m s _1 it remains to compute 
CDELT3Z. This is obtained by transforming CDELT3 from the 
topocentric to the barycentric frame with Eq. (8) expressed as 

A A C ~ V 

Av = Avp z , 

Vc 2 - v 2 

whence the frequency increment becomes 97.64775 kHz. 
Using Eq. (11) the corresponding wavelength increment is as 
recorded by CDELT3Z in Table 14. While SPECSYSZ is set to 
'BARYCENT', SSYS0BSZ records the fact that the observation 
was actually made from the ' TOPOCENT ' frame. 

The correction for the observatory's motion due to the 
Earth's rotation and orbital motion with respect to the barycen- 
ter amounted to 26. 108 km s" 1 in the direction of the source. As 
illustrated in Table 15, this velocity shift can be represented in 
the FITS header with the new keyword VEL0SYS. 

The barycentric frequency may be used as another alternate 
axis description. The reference frequency and frequency incre- 
ment for SPECSYSF = 'BARYCENT' have already been com- 
puted above. The relevant keywords have axis description F 
(for frequency) in Table 15. 

The frequency description may be translated into a wave- 
length description simply using Eqs. (10) and (11), shown in 
Table 15 using axis description W (for wavelength). 

The frequency description with respect to the BARYCENT 
may also be expressed simply as a radio velocity (also regu- 
larly sampled) using Eqs. (20) and (21). This is shown as axis 
description version R (for radio). 

A apparent radial velocity description (version V for veloc- 
ity) requires the use of Eq. (14) for the reference value and 
Eq. (15) for the increment. 

10.2. Real-world complications 

An assumption made throughout this paper, particularly in 
Sect. 7, is that data are provided with at least celestial coor- 
dinate information in addition to the spectral coordinates dis- 
cussed here. The normal method of providing such data is with 
the full set of coordinate keywords including WCSAXES with 
a value of at least 3 (for two celestial and one spectral axis) 
even if only one celestial position was observed. For spectra 
in tables, it would be possible to have simple columns for the 
celestial coordinates, labeled as such. But, as discussed at the 
end of Sect. 6.1.4, this simple case may still use the full WCS 
nomenclature in order to associate the coordinate columns with 
the data column. 

One important assumption made in Sect. 10.1 was that the 
topocentric velocity does not change significantly during the 
course of the observation. We now consider the consequences 
if this assumption is violated. In this we are interested in partic- 
ular with three-dimensional "data cubes" containing multiple 
samples on each of two spatial axes and one spectral axis. 

For some instruments, it may be possible to correct the data 
for each short exposure to a consistent WCS and then combine 
exposures to improve the sensitivity. For example, single-dish 
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radio telescopes observe a single direction at a time, allowing 
the spectral axis in each observed spectrum to be scaled and 
resampled before it is added into a three-dimensional image. 
For this type of instrument it is a relatively simple matter to 
correct for topocentric motion. 

However, other instruments observe multiple directions si- 
multaneously and also routinely combine multiple exposures 
to form an image. While the error over a single exposure may 
be negligible, images made by aperture synthesis radio tele- 
scopes, for example, are often based on multiple exposures 
taken months apart. These exposures are usually made with dif- 
ferent array configurations, which may change on a timescale 
of months, in order to sample the Fourier space adequately for 
imaging and deconvolution. Therefore, these data are affected 
not only by the rotation of the Earth but also by its motion 
around the Sun which usually is not negligible. 

As a specific example, the actual observation on which 
Sect. 10. 1 is based required an integration time of several hours. 
During this time the observing frequency was adjusted for the 
topocentric Doppler shift so that the barycentric frequency at 
the map center remained constant. This ensured that the spec- 
tral feature of interest, associated with a source near the map 
center, remained in the same spectral channel. However, this 
means that the primary spectral axis description, based as it is 
in the topocentric frame, is not strictly correct. It was used in 
this data as an aid to interferometric analysis, the resulting error 
being considered negligible, but strictly speaking a barycentric 
frequency reference frame, the F alternative WCS in the header, 
is the correct one. 

Even so, the barycentric frame is only strictly correct at 
the map center since the topocentric Doppler correction varies 
across the field-of-view; the correction is computed as the dot 
product of the observatory's velocity vector and a unit direction 
vector which, of course, varies across the field. This gives rise 
to a differential error that is greatest when the velocity vector 
is orthogonal to the direction to the source and least when it is 
parallel or anti-parallel. In other words, the differential correc- 
tion is greatest when the correction is zero, and least when the 
correction is greatest. 

For diurnal rotation the worst case differential error occurs 
for a telescope at the equator observing a source on the lo- 
cal meridian. For an angular offset, y, from the tracking cen- 
ter the maximum differential error is v sin y which amounts to 
8.7 ms -1 for y = 1°. The worst case change in the differen- 
tial error over the course of 12 hours occurs for an equatorial 
telescope observing a field at the north or south point on the 
horizon; it is 2usiny or 17.5ms- 1 for y = 1°. If data from 
exposures separated by several months are combined then the 
Earth's 30 km s _1 orbital velocity scales up the worst-case error 
to l.Okms 1 fory = 1°. Note that, in the example of Sect. 10.1, 
the barycentric correction is close to the maximum value so the 
differential correction is nearly at its minimum. 

The relevance of the SSYSOBSa keyword immediately be- 
comes apparent in this context. It records the reference frame 
in which the frequency has no spatial variation; in other frames 
the reference frequency and frequency increment are only 
strictly correct at the reference point. While SSYSOBSa will 
most often indicate the topocentric frame, it is possible to regrid 



(resample) the spectral data by applying a position-dependent 
shift-and-scale so that SSYSOBSa has some other value. 

Differential Doppler effects are often neglected for contem- 
porary observations made with small angular fields or modest 
spectral resolution. However, instruments now under develop- 
ment will be capable of significantly greater sensitivity over 
wider fields-of-view and/or spectral resolution and hence will 
need to consider the effects discussed here. 

11. Summary 

The new keywords required for spectral coordinate sys- 
tems are summarized in Table 13. SPECSYSa, SSYSOBSa, 
SSYSSRCa, and ZSOURCEa are used to specify velocity ref- 
erence frames; MJD-AVG, DATE-AVG, OBSGEO-X, OBSGEO-Y, 
0BSGE0-Z, VELOSYSa, and ZSOURCa enable conversion be- 
tween these standards of rest; and RESTFRQa, RESTFREQ, and 
RESTWAVa define the spectral line for which velocities are mea- 
sured. Variants of these keywords for use with tabular data are 
also defined. 

These new keywords and their allowed values, along with 
new values for the standard header keywords formalized by 
Paper I, and the associated algorithms and methods introduced 
here allow an accurate description of spectral coordinates in 
FITS images. Wavelength, frequency, apparent radial velocity 
and oft-used quantities that are linearly related to them, such as 
redshift, energy, and radio velocity, may now be expressed as 
functions of pixel coordinate along an axis regularly sampled 
in wavelength, frequency, or apparent radial velocity. 

A non-linear algorithm is also provided for spectral dis- 
perses commonly used at optical wavelengths: gratings, 
prisms, and the combination of the two, grisms. Although de- 
veloped for ideal dispersers, it was shown to apply quite well 
to real-world dispersers with suitable fine adjustment of the in- 
strumental parameters. 

The multi-dimensional -TAB table lookup algorithm devel- 
oped for wavelength calibration will be useful for much more 
than just spectral axes. Several others of the methods intro- 
duced here are also completely general; logarithmic (-LOG) 
coordinates, and the CNAME ia keyword. The mathematical for- 
malism described in Sect. 3.4 for constructing one-dimensional 
non-linear coordinate systems also has wide validity. 

Appendix A: Relativistic space velocities 

The discussion of Sect. 2, particularly Fig. 1, indicates the im- 
portance of transverse, as well as radial velocities, in spec- 
troscopy at significant velocities. To be concrete, consider a 
relativistic jet emanating from a distant galaxy. In principle, the 
galaxy's systemic, cosmological redshift can be measured sep- 
arately and used to correct the jet's observed redshift thereby 
providing its kinematic redshift (i.e. associated with a true ve- 
locity) in the reference frame of the galaxy. 

Clearly knowledge of the velocity is fundamental in study- 
ing jet kinematics and dynamics. However, it can only be com- 
puted from the kinematic redshift if the jet's orientation angle 
is known. Note that the equations involving velocity in Table 3 
are actually only correct if the transverse velocity is zero. For 
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Table A.l. Velocity equations using orientation and total velocity. 



V = VQ 



y^ 2 - j| 

c + u s sin # ' 



dv u s + c sin 

T~ = ~ cy « , 

df s ( c + v s sin £>) 2 y c 2 - « 2 



c + y s sin 6 

d/1 u s + c sin 

d^ = ^V-f 2 ) 3 ' 2 ' 
-y 2 sin ± iff 
X 2 + <P 2 

where x = — = — , 

v /i 



and iff = ^jl - x 2 cos2 

du s 1 du s 

dv v dx ' 

dv s -A dv s 

dA ~ ~dx' 

dv s 2^sin6»±(l + ^ 2 sin 2 9) 

where ¥ = "^ ^¥) 2 



(A.4) 
(A.5) 

(A.6) 
(A.7) 
(A.8) 

(A.9) 
(A. 10) 
(A.ll) 
(A. 12) 
(A.13) 



example, Eq. (12) is just a special case of Eq. (2). This is why 
the velocity is always referred to in the paper as being "appar- 
ent." 

There are instances where the orientation angle may be in- 
ferred by geometry (e.g. by the observed tilt of an accretion 
disk) or by modelling. Defining the orientation angle as 8 we 
have 



i> r = v s sin ft, 
u t = v s cos 8, 



(A.l) 
(A.2) 

(A.3) 



where v t is the radial velocity, v t is the transverse velocity, and 
v s > is the total space velocity. The orientation angle 8 is 
defined so that 8 = -90° is towards the observer, 8 — is trans- 
verse, and 8 = +90° is away from the observer. It is beyond 
the scope of this paper to describe how 8 should be determined 
or what exactly it means for particular relativistic observers, 
other than that it resolves the kinematic velocity into radial and 
tangential components. With these definitions, the equations of 
Table 3 are given in Table A.l. Note that for^ > 1 Eqs. (A.8) 
and (A.13) have two valid solutions for 8 < but none other- 
wise. For^f < 1 they have only one valid solution for any value 
of 8. Invalid solutions, i.e. those with v s < 0, are the negative of 
the valid solution for -8. 

There are other combinations that are of interest, including 



sin# 



vo 



vv s 



c 
c 



It would perhaps be more instructive to express these sorts 
of equations in terms of the actual radial and transverse veloc- 
ities. This is relatively straightforward, but the equations be- 
come even more messy and so we omit them here. 

It has been proposed that we introduce a single keyword 
to express 8 so that we may express velocities in terms of a 
true space velocity u s rather than an apparent velocity u as done 
throughout this paper. However, astrophysical jets are thought 
to be conical in form, to interact along their boundaries with the 
external medium, and to change direction, even assuming he- 
lical shapes. Thus, there is no single angle 8 that can normally 
be used to describe all directions in an image and one should 
probably not even assume that 8 is independent of redshift in 
any one direction. For those cases in which a single angle may 
be appropriate, such as a spectrum of a spatially limited region 
in a relativistic jet, we reserve the keyword 



VELANGLa 



(floating-valued) 



with default value +90° to express the angle 8 in degrees. All 
of this discussion has neglected non-velocity redshifts such as 
those due to the gravity of the black hole thought to be at the 
base of the astrophysical jet. Therefore, we choose to leave the 
full expression of internal velocities within a celestial object 
to those observers who are armed with both extensive obser- 
vations and a detailed model. Apparent velocity is a suitable, 
semi-physical concept appropriate to a general FITS standard. 
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